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VERTICAL DISTRIBUTION OF VIBRATIONAL ENERGY OF 
MOLECULAR NITROGEN IN A STABLE AURORAL RED ARC 
AND ITS EFFECT ON IONOSPHERIC ELECTRON DENSITIES 

by 

George P. Newton 
ABSTRACT 

Previous solutions of the problem of the distribution of vibrationally ex- 
cited molecular nitrogen in the thermosphere have either assumed a Boltzmann 
distribution and considered diffusion as one of the loss processes or solved for 
the energy level populations and neglected diffusion. This work combines both 
of the previous approaches by solving the time dependent continuity equations, 
including the diffusion process, for the first six energy levels of molecular ni- 
trogen for conditions in the thermosphere corresponding to a stable auroral red 
arc. The primary source of molecular nitrogen excitation was subexcitation, 
inelastic, collisions between thermal electrons and molecular nitrogen. The 
reaction rates for this process were calculated from published cross section 
calculations. The loss processes for vibrational energy were electron and 
atomic oxygen quenching and vibrational energy exchange. The coupled sets of 
non-linear, partial differential equations were solved numerically by employing 
finite difference equations. The results show that molecular nitrogen in excited 
vibrationally to the degree that the rate constant for the rate limiting ionospheric 
loss process 0 + + N 2 - NO + + N is increased by as much as a factor of eight at 
F 2 region altitudes. It was found that deviations from a Boltzmann distribution 
caused the reaction rate constant to be as much as 1.7 times the reaction rate 


constant calculated for the energetically equivalent Boltzman distributed densities. 



The time constant for loss of ionospheric electrons was calculated using the en- 
hanced reaction rates. It was found that, for altitudes where the normal loss 
time constant is less than one day, the loss time constant was decreased in the 
arc by factors of 1.6 to 5 depending, on altitude. It is concluded that the ob- 
served decrease of electron density in the F 2 region in stable auroral red arcs 
can be explained by enhanced reaction rates for ion-atom interchange between 
0 + and N 2 caused by vibrational excitation of molecular nitrogen through 
electronic collisions. 
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CHAPTER I. INTRODUCTION 


1.1 Description of Stable Auroral Red Arcs 

Stable Auroral Red Arcs (SAR-Arcs) are photometrically observed enhance- 

o o 

ments of air-glow at 6300 A and 6364 A wavelengths occurring during geomag- 
netically active times. They were first observed by Barbier (1) in 1957. The 
emissions are from the forbidden red lines of atomic oxygen OI ( 1 D- 3 P). Peak 
intensities of emissions are located at altitudes near 400 km. and vary from a 
few hundred rayleighs* to tens of kilorayleighs. The arcs generally occur near 
midlatitudes and are distinct in location from usual auroral phenomena. They 
extend over spatial dimensions of approximately 1000 km. horizontally in a 
meridional direction and globally in a longitudinal direction. The arcs are 
magnetically controlled as evidenced by their alignment in geomagnetic, rather 
than geographic latitude, and by their simultaneous occurrence at magnetically 
conjugate points. Once the SAR-arc is formed it generally persists for ten hours 
or longer. The frequency of observation appears to follow the solar cycle. Re- 
views of the general properties of SAR-arcs and citations to the literature are 
contained in Reference (2) and more recently Reference (3). 

1.2 Cause of SAR-Arcs 

The absence of the atomic oxygen green line OI ( 1 S- 1 D) in the air glow en- 
hancements during SAR-arcs, implies that a low energy source must exist for 
the excitation of the OI( x D) state. Several mechanisms have been proposed for 
selectively exciting X D level and their implications studied (see (3) for summary) 

*a rayleigh is 10 6 photons cm -2 sec -1 
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but one proposal has recently gained prominence. In this proposal (4) energy 
originating in the solar wind is transferred through the ring current in the mag- 
netosphere to the ambient plasma, from which it is conducted down magnetic field 
lines into the ionosphere where it heats ambient electrons. These heated electrons 
collisionally excite the atomic oxygen *D level, but lack sufficient energy to excite 
the *S level. This proposal was utilized together with in- situ measurements of 
electron temperatures and densities by Roble et al. (5), (6) who found it to be con- 
sistent with all observations. Other proposals have been excluded by in-situ 
measurements (see (3) for discussion of the proposals and (7), (8) for observations). 

1.3 Ionospheric Behavior in a SAR-ARC 

From ground based (see (3) for review) and more importantly from satellite 
obtained data (9), (6), (7), (8), the following behavior of the ionosphere in the vi- 
cinity of SAR-arc has been determined. 

At the satellites where the electron temperatures, and in some instances the 
ion temperatures, were measured a general enhancement of the charged particle 
temperatures was observed on the geomagnetic field lines intersecting the SAR- 
arcs. The electron densities at these satellite attitudes (usually greater than 
1000 km.) showed for different SAR-arc occurrences both increases, decreases 
and no change in the electron densities, apparently depending on the altitude of 
measurement. However, corresponding topside sounder data always showed 
electron density depressions in the lower F 2 regions (near 300-800 km altitude) 
in the corresponding region of enhanced electron temperature. The electron 
density scale heights often are less in the red arc region, reflecting changes in 
the ion composition (7), (6). The more complex electron density behavior at the 
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higher altitudes (1000 km and greater) may in part be explained by a combination 
of the ion composition changes and the increased plasma temperature. A typical 
observation of the electron density and temperature in the vicinity of a SAR-arc 
is shown in Figure 1-1. 

1.4 Objective of this work 

One of the primary loss processes for ionization in the F 2 region, where 
atomic oxygen ions are the dominate ion species, is 

0 + + N 2 i NO + + N (1-1) 

followed by 


NO + + e i N + 0 (1-2) 

where k and a are the ion-atom interchange and dissociative recombination re- 
action rates, respectively. Typically a is about four orders of magnitude greater 
than k, (11), so that the rate determining process for ionization loss is reaction 
(1-1). However, k in reaction (1-1) is known to be a strong function of the vi- 
brational temperature of molecular nitrogen (12), (13). This is why knowledge 
of the energy distribution of the excited levels plays a crucial role in the under- 
standing of the electron loss rates in the SAR-arc. Thus, as suggested earlier, 
(12), the loss of ionospheric ionization may be considerably enhanced by the 

presence of vibrationally excited nitrogen. The source of the vibrationally ex- 

2 

cited nitrogen molecules is sub- excitation, inelastic, collisions between the 
2 

Subexcitation means that the final state of the N 2 is neither electronically ex- 
cited nor ionized. 
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Figure 1-1. Observed ionospheric electron densities and observed and calculated electron temperatures for a 500 R SAR-arc 
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Figure 1-2. Model neutral atmosphere and ionospheric electron densities 
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Figure 1-3. Neutral atmosphere kinetic and electron temperature models 
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ambient molecular nitrogen and the thermal ionospheric electrons, the same 
mechanism that excites the OI^D). The cross section for this type of process 
is important to this study (Appendix I). 

The objective of this work is to investigate the distribution of vibrationally 
excited molecular nitrogen under conditions corresponding to a SAR-ARC to 
determine if the presence of vibrationally excited nitrogen can account for the 
observed decrease in the electron density in the F 2 region. For this study we 
used the same electron density profiles that Roble (10) used to calculate the 
temperature profiles for SAR-arcs, and we used his calculated electron tem- 
peratures. The electron density profile and model atmosphere used for this 
work are shown in Figure 1-2. The neutral atmospheric model used for this 
work is the 1965 Jacchia model (14) with an exospheric temperature of 1500°K. 
We did not include minor constituents. The neutral atmosphere temperature 
profile is shown in Figure 1-3, together with several of Robel’s, (10), electron 
temperature profiles corresponding to different SAR-ARC intensities. 



CHAPTER H. FORMULATION OF THE PROBLEM 


2.1 Previous Calculations of the Distribution of Vibrationally 

Excited Nitrogen 

Several authors have suggested that vibrationally excited nitrogen may in- 
fluence the excitation of atmospheric radiations, chemical reactions, and the 
ionospheric temperature (15), (16), (17), (18), (12), (19). Walker et. al. (20) were 
the first to calculate the altitudinal distribution and diurnal variation of the vi- 
brating molecules assuming a Boltzmann distribution. The diurnal calculations 
assumed an unperturbed atmosphere and ionosphere. They considered diffusion 
of the vibrationally excited nitrogen vertically and limited their calculation to alti- 
tudes between 100 km and 300 km altitude. They found that vibrational energy 
exchange with C0 2 was the dominant loss process for vibrational quanta below 
125 km altitude. The principal sources of vibrational quanta in the lower thermo- 
sphere (100 km and 600 km) were found to be: (1) the reaction between atomic 
nitrogen and nitric oxide; (2) the collisions of photoelectrons and nitrogen mole- 
cules; and (3) possibly the quenching of metastable oxygen atoms, OI ( X D), by 
nitrogen, although the vibrational yield of the last process is unknown. 

Recently, McNeal, et. al. (21) measured the reaction rate for vibrational- 
translation energy exchange between atomic oxygen and vibrationally excited 
molecular nitrogen at 300°K and found the rate to be considerably larger than 
extrapolated values (which were used by Walker et. al. (20).) obtained from ex- 
periments at higher temperatures. The quenching of vibrationally excited mole- 
cular nitrogen by atomic oxygen, which was ignored in the calculations of Walker 
et. al. (20), has been included in the calculations of Breig, et. al (22) who used 
the recently determined reaction rate. They essentially repeated the calculations 
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of Walker, et. al. (19) and found that inclusion of the loss of vibrational excitation 
by nitrogen molecules to atomic oxygen significantly decreased (by factors of 2 to 
4) the vibrational temperatures in the lower thermosphere. 

Schunk and Hays (23) solved the continuity equations for the first eleven energy 
levels of molecular nitrogen for altitudes above 150 km. under conditions corre- 
sponding to an aurora. Diffusion processes were not included in their calculations. 
The source of vibrational excitation was non-thermal electrons and the only loss 
process included for each level was vibrational energy exchange between nitrogen 
molecules. They found that when the non-thermal electron source was acting, the 
N 2 vibrational distribution was non- Boltzmann and that the deviation from the 
Boltzmann distribution increased with increasing altitude. They found that the 
loss rate for 0 + due to reaction (1-1) was increased by a factor of 1.5 when al- 
lowance was made for the non- Boltzmann distribution. The non-thermal electron 
source was allowed to operate for 15 min. and they found that the enhanced loss 
rate persisted for 1000-2000 sec. after the auroral bombardment commenced. 

Walker (24) discussed the result of Schunk and Hays (23) and included in his 
discussion the effects of diffusion. He starts with the steady state continuity equa- 
tion for each level and ignores quenching processes for the vibrationally excited 
nitrogen molecules. Under these conditions, the source of vibrational excitation 
and diffusion are the only terms which can cause departures from a Boltzmann 
distribution. By comparing the production term for a given vibrational level with 
the production at that level due to vibrational exchange, he concluded that de- 
partures for vibrational levels containing a significant fraction of the nitrogen 
molecules will not exist at altitudes below 300 km. By comparing the diffusive 



flux term with the loss due to vibrational energy exchange he concludes that dif- 
fusion is less important than vibrational energy exchange at altitudes below 370 
km. These conclusions were based on conditions of an undisturbed atmosphere 
and source terms similar to those used by Walker, et al. (20). 

2.2 Differences between this and previous works. 

The source of vibrational excitation for this work is the thermal electrons 
in SAR-arcs. This source extends to higher altitudes than the sources used by 
Walker, et. al. (20), Breig et. al., (22) and Schunk and Hays (23). As a result, the 
molecular excitation occurs in a region of the thermosphere where vibrational 
energy exchange is slow in comparison with the diffusion process. Thus, be- 
cause the source is non- Boltzmann in its energy distribution and because the 
excited molecules diffuse to lower altitudes where vibrational energy exchange 
is the dominant process, the assumptions of a Boltzmann distribution and negli- 
gible diffusion loss are no longer tenable. This work differs from that of Walker, 
et. al., (20) and Breig, et. al. (22), by not assuming a Boltzmann distribution, and 
from that of Schunk and Hays (23) by including diffusion of the excited molecules. 

2.3 The equations of continuity. 

The continuity equation for each energy level of vibrationally excited 
molecular nitrogen, N* , governs the level population and its altitudinal distribu- 
tion. For the ith level the equation is 
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where n. is the number density of molecular nitrogen in energy level i 
v. is the vector velocity of transport of n. 

Q L is the production rate of m 
L. is the loss rate of n . . 

X 1 

In principal the production and loss rate terms contain all sources and sinks for 
the ith level. However, we are interested in solutions of (2-1) corresponding to 
SAR-arc conditions where the electron temperature is unusually large. In addition 
we are interested in a process which, to first order, can be considered as an ad- 
ditional source to those ordinarily present in the thermosphere. 

Figure 2-1 shows the rates of the dominant processes for equation (2-1). 

Also included in Figure 2-1, for comparison purposes, are the total vibrational 
quanta production rates between 120 and 300 km. calculated from the work of 
Walker et. al. (20) and Kummuler and Bortner (25) and shown as squares and 
circles, respectively. Shown as triangles in Figure 2-1 are the C0 2 quenching 
rates for vibrationally excited nitrogen at 120 km. taken from references (20), 

(22) and (25). 

The atomic oxygen quenching curve shown in Figure 2-1 was calculated using 
the reaction rate of 3.5 x 10 -ls cm 3 sec -1 , (21), and the atomic oxygen number 
density shown in Figure 1-2. The curve of atomic oxygen production of N * (trans- 
lational-vibrational energy enchange process) was calculated from the quenching 
curve using detailed balance considerations. These curves correspond to 
1-0 and 0-1 transitions of N 2 respectively. 

The vibrational energy exchange curve was calculated for the 0 - 1, 1 - 0 
transitions using the exact resonance case of the Schwartz-Slawsky-Herzfeld, 
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(SSH), theory, (27) page 330, and the molecular nitrogen number density profile of 
Figure 1-2. 

The rate for diffusion was calculated from the diffusion coefficient, D, and the 
density scale height, H, for molecular nitrogen by 

rate = -P_ (2-2) 

H 2 

This rate may be regarded as the inverse of the time for a nitrogen molecule to 
diffuse to an altitude where there is an appreciably greater density (by a factor of 
e) of molecular nitrogen. 

The electronic production curve was calculated using the reaction rates for 
inelastic collisions between thermal electrons and molecular nitrogen discussed 
in Appendix I and the electron and molecular nitrogen densities and electron tem- 
peratures in Figures 1-2 and 1-3, respectively. The electron temperature used 
is that for a 17.7 KR SAR-arc. For the calculation it was assumed that the vi- 
brational temperature was zero. 

From Figure 2-1 it is seen that the source of vibrational energy from 
thermal electrons, under the SAR-ARC condition, is larger over most of the alti- 
tude range of interest than the source of vibrational energy from ordinary 
ionospheric processes (20), (25). Because of the uncertainty in the yield of vibra- 
tional energy from 0 ^D) quenching by molecular nitrogen (20), (25), and because 
the O (*D) density during SAR-arcs is enhanced and its spatial distribution is 
complex, we have ignored this potential source term. Thus, the source term, 

Q. , in equation (2-1) will only include production of vibrationally excited nitrogen 
by collisions between nitrogen molecules and thermal electrons (and only 
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coincidentally by atomic oxygen) and vibrational energy exchange between nitro- 
gen molecules. These approximations should underestimate the source terms and 
hence the calculated temperatures would tend to be a lower bound for this reason. 

The dominant loss processes for vibrationally excited molecular nitrogen at 
higher altitudes is electron quenching at altitudes above 400 km. (20, 25), and 
diffusion above 320 km. (see Figure 2-1). At low altitudes the dominant loss 
processes are quenching by atomic oxygen for altitudes between 120 and 300 km 
(Walker et. al., (20), found that quenching by other minor constituents in the lower 
thermosphere through translational-vibration energy exchange was negligible) and 
vibrational energy exchange with C0 2 (20), (22), at 120 km. Thus, the loss term, 
L., in equation (2-1) will include loss by electron and atomic oxygen quenching and 
vibrational energy exchange with molecular nitrogen. The vibrational energy loss 
by molecular nitrogen to carbon dioxide will be included through a lower boundary 
condition on the vibrational temperature for the following reasons. The abundance 
of carbon dioxide at 120 km is somewhat uncertain and depends on the altitude of 
the turbopause. The turbopause of the atmosphere is below 120 km altitude. Thus, 
in the altitude region of interest, carbon dioxide could be expected to approach dif- 
fusive equilibrium. Because of its relatively large mass, the density of carbon 
dioxide should decrease rapidly above 120 km. and increase rapidly below 120 
km. altitude, thus creating a loss process significant only at the lower boundary. 
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where a prime on the summation excludes j = i, and 

**i , j j i-i , j +1 * s exc h an S e transition probability for a nitrogen molecule in 
level i - 1 going to level i while the collisional partner molecule 
goes from level j + 1 to level j. 

za . is the collision frequency of a nitrogen molecule in level j with a 
molecule in level i. 

C. . is the probability for collisional transitions from level j to level i 
where the collision is with atomic oxygen atoms. 

A. . is the reaction rate, from Appendix I, for the transition from level 
j to level i due to collisions between thermal electrons and molecular 


nitrogen molecules in level j . 
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v k [ 0 ] is the collision frequency of a molecule in level i with an atomic 
oxygen atom. 

2.4 The coefficients in the equations of continuity. 

The cross sections for the subexcitation, inealstic, scattering of electrons by 
molecular nitrogen (28), (29), (30) are only available for the first five energy levels. 
These cross-section calculations and the calculation of reaction rates from the 
cross-sections are discussed in Appendix 1. 

We will assume that molecular nitrogen can be represented by a harmonic os- 
cillator. Thus, only transitions between adjacent levels will be considered. This 
assumption is appropriate for the moderate vibrational temperatures of this prob- 
lem where only the lower lying energy levels are important and where anharmonic 
effects are small (31), (32), (33). It is also implicit to this calculation that the first 
five excited levels of molecular nitrogen will be considered as the totality of excited 
levels. 

The above assumption and the theory for exact resonance vibrational-vibrational 
energy exchange of SSH, (27), allow us to write the vibrational-vibrational energy 
exchange probability in terms of the 0 - 1, 1 - 0 transitions of molecular nitrogen 
as 


^i, jsi— l.j+I ~ + ^*1, 0:0,1 


From detailed balance considerations 


P = P 

i-l , j +1 : i , j 


(2-5) 


i , j : i~l . j + 1 


(2-6) 
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Setting the steric factor in the SSH theory to three, and using the Lennard-Jones 
potential parameters for N 2 from table 66-1 of (27), and table 66-5 of (27) for the 
characteristic % of the theory, we have 

p i.j,i-i, j+ i = 2 .6°°322 x 10-6 i(j + 1) T e 91,s/T (2-7) 

where T is the kinetic gas temperature. 

It is assumed that 


c i+i.i = + !) C i.o (2-8) 

for the translational-vibrational energy exchange probability, the reason will be- 
come apparent later. Thus, (2-8) is consistent with SSH theory, 3 but C x 0 will not 
be determined from this theory (see reference (21)). Detailed balance gives 

C i.i-i =c i-i,i e ~ 8/T (2-9) 

where 6 is the characteristic temperature of N 2 (d = 3380°K). 

We assume that the vibrational level in which the nitrogen molecule resides 
does not effect its collision cross section for both molecules and atoms. Thus we 
write 


v . . = v n. 

j i i 


v i.[o] = v ' n[0]. 

Equations (2-3) and (2-4) become respectively 

Qi = Z] "i + P l. 0:0.1 V YL i(j + 1} n i-l "j+1 + P l. 0:0,1 V 

j=° j=0 

YL j (i + 1} n i+i n j-i + C 1.0 v ' n [0] [i n i-i + (i + 1) e 5/T n. +1 ] 


j = i 


(2-10) 

( 2 - 11 ) 


(2-12) 


The measured reaction rate for C Q>1 is several orders of magnitude larger than 
would be calculated using SSH theory. The reason for the discrepancy may be that 
chemical interactions are responsible for the anomalous quenching, whereas SSH 
theory attributes the interaction as being due to short range repulsive potentials. 
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j i 


n. + 


1 , 0 : 0,1 


V 


21 * 


i(j +1) n. n + P 


1 , 0 : 0,1 


j=o 


y J (i + 1) rij n. + C l 0 v' n [0] [i e 0/T + (i + 1)] n. 
j“i 


(2-13) 


The quantity v in equations (2-10), (2-12) and (2-13) is calculated from equation 
(36-3) of reference (27) 


V 


1.271 p 
^ 2 n 


(2-14) 


where 

p is the gas pressure 
N 2 is the viscosity for molecular nitrogen 
n is the total number density of molecular nitrogen 
Using the ideal gas law, we have 


V - 


1.271 kT 


(2-15) 


where k is Boltzmann's constant. Using (8.2-18) of reference (34), and average 
values of the appropriate parameters for temperatures of 355°K and 1500°K, we 
obtain for v 


v - 1.482094 x 10” 11 /T (2-16 

The quantity Cj 0 v‘ occurring in equations (2-12) and (2-13) is calculated 
from C 01 v • by equation (2-9), where C Q1 v' is considered to be a 'lumped" con- 
stant equal to the reaction rate constant reported by McNeal, et. al., (21) and is 
considered constant with altitude. 

It will be assumed that the neutral atmosphere is static, and uniform hori- 
zontally in all of its parameters. The transport velocity, v*. , is then the diffusion 
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velocity for molecular nitrogen through the atmosphere, and is given by (35), page 
244, 

Vj = - D ["— (Vp. - p, g) + k T 1vt"1 (2-17) 

1 p. 1 x T 

_ 1 -J 

where 

D is the binary diffusion coefficient for molecular nitrogen, assv med 
to be independent of vibrational excitation, 
p. is the pressure of the molecular nitrogen in energy level i. 
p is the mass density of molecular nitrogen 
k T is the thermal diffusion ratio 
g is the acceleration of gravity. 

Using the ideal gas law in (2-17) we have 

v i= -D iVn. + (k T + l)IvT + £|k (2-18) 

where k is a unit vector in the vertical direction 


m is the mass of molecular nitrogen. 
Rewriting (2-18) we have 


where 



I = f k + 1)1 — +HI 

H ( T ) T dz kT’ 


(2-19) 

( 2 - 20 ) 


and z is the vertical coordinate, positive upward and H is the density scale height. 
Thermal diffusion for molecular nitrogen is small, and therefore k T will be set to 
zero. The flux term in equation (2-1) becomes 



( 2 - 21 ) 
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It will be assumed that the magnetic field line is vertical, so that the 
SAR-arc is symmetrical about the vertical. Under this assumption the hori- 
zontal flux at the arc center vanishes. If, in addition, the divergence of the 
horizontal flux is ignored, equation (2-1) can be reduced to one spatial dimen- 
sion corresponding to conditions at the arc center. 

Under the above assumptions and using (2-21), (2-12) and (2-13), equation 
(2-1) becomes 


d n 


‘+ #- (F 


d t dz 


i),=£ {A i j n j ~ A j i n i> + P l. 0:0,1 Vi Yl ( J +1) [n i-l n j +1 

i - 0 ; = 0 


-nm.] + Pj, o:o,i v + X ) YL ] [n i + i n j -i " n i n j 


(2-22) 


C 1>0 v' n [0] [i n._ 1 + (i + l)e e/T n i + 


- (i e 6/T + i + 1) n £ ] 


2.5 Boundary conditions. 

Two sets of boundary conditions will be used for equation (2-22). 
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For the first set we have 

Initial condition: n i (z, 0) = 0 for all i > 0 

= n for i = 0 


Lower boundary: n. (120, t) = 0 for all i > 0 

= n (120) for i = 0 


Upper Boundary: 



n £ (®, t) . 

H 


for all i . 


The initial and lower boundary conditions correspond to all molecular nitrogen 
being in the ground state. This would correspond to total quenching by C0 2 at the 
lower boundary, and no excitation initially. The upper boundary condition insures 
zero flux, for all excited levels, at the top of the atmosphere. 

For the second set of boundary conditions, 


Initial Condition n.(z, 0) = 2 e -i 0/355 for all i 

1 1 - e ~ 0/355 

Lower Boundary n . (120, t) = 2 e -i 5/355 for all A 

1 1 _ e - 0/ 3S5 

Bn. n (co, t) 

Upper boundary = for all i 

r=co H 

The initial and lower boundary conditions in this case correspond to assuming the 
vibrationally excited nitrogen to be Boltzmann distributed corresponding to a vi- 
brational temperature of 355°K, the kinetic temperature at the lower boundary. 


2.6 Simplification of equations and technique of solution. 

In Appendix II the set of six, coupled, non-linear, partial differential equations 
(2-22) are replaced by seven equations in which the essential coupling is more 
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readily observed. This simplification is achieved by expressing the solutions of 
equation (2-22) as a sum of a Boltzmann distributed part and a part representing 
a deviation from a Boltzmann distribution, that is 

n. = n! + e. (A2-1) 


It is then required that the Boltzmann distributed part, n ' , possess the total vi- 
brational energy and number density of molecular nitrogen for the problem (see 
equations (A2-2) and (A2-3)). The deviation part, e i , is then a measure of the de- 
parture, at each altitude, of the solution from the energetically equivalent Boltz- 
mann solution part. The differential equation to be solved for the Boltzmann part 
of the solution is 


B0 l_3_DnBi/, p gV 
Bt n 3z Bz g z 2 


= H H Ijti (fnj 4] 

j * 0 i = 0 


(A2-30) 


+ Cj ^ 0 v* n tO] [i/< + 1 - e 9/T i//] 


where i p is related to the partition function, Z p , for the vibrationally excited 
nitrogen by 


0 = (Z p - 1) 


(A2-27) 


The differential equation for the deviation part of the solution is 
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i + 1' 

1 

Bqp. 

U« + l 
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1 

J" 

0 + 1. 
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"1 Bi p B<// 

/i 

i 
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0 + 1 
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10 

"0 
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(1 

i 

+ il 

B 2 0 
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0 + 1J 

tBz // D 

10 

“0 

+ ij 

’ Bz 2 

\ ^ z / 


* T Aji (1 - s i|) + p ., 0,0.1 "" « 2i *')**« 
(^ + 1)J ^o J n 


+ Cj f 0 v' n [0] ( i e 6/T + i + 1 ) 


cp i = 2I{ A ij (i+<p i )e " (i ’ 1> ' /Tv } (i - 8 ii ) 

j = 0 ^ ■* 


j_ L±l \ IhL ]tL fi_ j + i \ f3D r _i_ 

0 + ly Bt Bz \i fj \p + lj l_Bz + \ H 

/* mxitl* rLii-M-nILMi 

\0 0+lJ ^zj Z— j J1 lJ \0 0 + 1/ 3 z 2 


- D 


' M 2 f i . i + 1 

, BzJ 


+ ~-~' 2 } +P 1. 0:0.1 Vn{ ( 2i + 1)0+ 0 

Y (0 + 1) J 


+ C 10 v'n[0] (i e 6/T + i + 1) 


+ P l. 0:0.1 pn (< 2i + 1) 0 + i> 


+ C 1>0 v'n[0] [i e^v + (i + 1) e^-^/Tv] + Vi ^ [ P i. 0 : 0 .i ™ i V 
+ C l(0 v'n[0]i] +cP i+1 e‘ 5/T v [P lt0s0<1 vn(i +1) (0 + 1) 


+ Cj Q v' n [0] e 0/T (i + 1)J 


(A2-46) 
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where 


e. 



The boundary conditions are 
set Is 


(A2-39) 


Initial condition 

0(2 

:, 0) = 0 


Lower boundary 

0(1 20, t) = 0 


Upper boundary 

Bi p 
Bz 

N 

II 

8 

1! 

O 

(A2-31) 

set 2: 




Initial condition 

0(2 

:, 0) = 7.329241 x 10" s 


Lower boundary 

0(120, t) = 7.329241 x 10“ 5 


Upper boundary 

Bi// 

Bz 

o 

II 

8 

II 

N 

(A2-32) 

set 1 and set 2 




Initial condition 

<Pi( 

z, 0) = 0 foralli 


Lower boundary 

<p.(120, t) = 0 for all i 


Upper boundary 

Bcp. 

Bz 

= 0 foralli 

i=0 

(A2-47) 


The set 1 boundary conditions on 0 correspond to a vibrational temperature 
initially zero, with the vibrational temperature at the lower boundary remaining 
zero, presumably due to quenching. The set 2 boundary conditions for 0 correspond 
to the initial vibrational temperature of 355 °K (the kinetic temperature at the lower 
boundary) and the vibrational temperature at the lower boundary remaining at 
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355 °K. In both cases the upper boundary condition on p corresponds to a zero flux 
of quanta through the top of the atmosphere. 

The bo undar y conditions on the cp. ’s correspond to zero deviations from a 
Boltzmann distribution initially, with zero deviation being maintained at the lower 
boundary (presumably due to the rapid vibration-vibration energy exchange among 
N 2 at this altitude). 

Calc ulati on shows that the source term for vibrational quanta due to electron 
collisions in equation (A2-30) is only mildly non-linear (less than 30%) for the elec- 
tron an d vibrational temperatures of interest here. It can also be argued, that for 
the relatively weak source of vibrational excitation for this problem, the departures 
from a Boltz mann solution will be such that the term e . /n will be small relative to 
1. Thus by dropping the e j /n term in (A2-30), the P and cp* equations are un- 
coupled, and a first approximation to p may be obtained by solving the equation 


"bp 1 BDn bp 
Bt n bz bz 


Bz 2 


" 1 

( + VI 

_sT+l 

[pnll 


(i - i) A u 


+ C li0 v'n[0] &/> + ! -e 0 ' r p] 


(2-23) 


When p has been obtained, the vibrational temperature may be calculated by 


T v (z, t) = 


e 


In 


0(z, t) + f 
p(z, t) _ 


(2-24) 


and the coefficients of equations (A2-46) are known functions of altitude and time. 
A first approximation to the cp. ’s can then be obtained by solving equation (A2-46). 
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If necessary, a second approximation to 0 may be obtained from equation 
(A2-30) by using the first approximation of the cp. ’s to determine e /n in the source 
term of equation (A2-30) and resolving equation (A2-30). A second approximation 
to the qp ’s can then be obtained by using the second approximation for v p in equations 
(A2-46). Higher approximations may in principle be found by repeating this process. 

2.7 Finite difference equations and solutions. 

In Appendix HI it is shown that the finite difference equations for both equations 
(A2-30) and (A2-46) take the general forms 


and 


- A 


i K\ + B 4 C - m 


(A3- 7) 


-*A 


^ +1 cp m+1 + ^B. 
^1 1 


'fyn+l 


1^1 1 


(A3-17) 


where 1 refers to an altitude grid point, with 4=1 the bottom boundary and 4 = L the 
top boundary, m refers to a time grid point and i refers to the energy level being 
solved for. The boundary conditions become for p 


set 1 


set 2 


m = 1 0* = 0 

(A3-8) 

4 = 1 </»” = 0 


= 7.329241 x 10" s 

Jjtn 0 1 

“ 4 


(A3-9) 


4 = L 


'A” - K-i = 0 


'PI - 'PI-! = 0 




and for the ^ 's 
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m = 1 

** 

1 = 0 

for 

all i 


i = i 


? = 0 

1 

for 

all i 

(A3-18) 

l = L 

L 0' 

■ - L ~ 1 d) m = o 
1 ^1 

for 

all i 



The two equations (A3-7) and (A3-17) are of the same form. We follow Richtmyer, 
(36), pages 103-104, for the solutions of (A3-7) and (A3-17). Only the solution of 
(A3-7) is presented here, because the extension to each of the equations in (A3-17) 
is clear. 

It is seen that the right hand side of (A3 -7) is a known quantity. Equation (A3- 7) 
is a one parameter family, because any value of determines a solution. Let 

(2_25) 

which is also a one parameter family. We seek to find and such that the two 
families of solutions (2-25) and (A3-7) are the same. If (2-25) is to be true for any 
member of the family, the boundary conditions (A3-8) and (A3-9) require 
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n 


m + 1 _ 




'P 


m + 1 


°l + H F -t- , 


B^? — . ^+1 Bp— A^Ej 


(2-29) 


V " A «-i 


Comparing equation (2-29) with (2-25) we have 



_h 

~ h H- 1 



D i + h ^c-i 

E ^-i 


(2-30) 


(2-31) 


Using (2-26) and (2-27) or (2-28) we can calculate the E^'s and F^’s from 1 = 1 to 
1 = L, that is, from the lower boundary to the upper boundary, by using (2-30) and 
(2-31), respectively. 

When 1 = L - 1, equation (2-25) gives 




ii m + l 
L 


+ F, 


L- 1 


(2-32) 


Using the upper boundary condition of (A3-8) or (A3-9) in (2-32) we have 


0 


ra + 1 
L- 1 


L- 1 


>-Vi 


(2-33) 


We now repeat the process, for ^ = L - 2, L - 3 1, to solve for the p m+ 1 's, 

that is, we solve for the '/ ,m+ x, s from the upper boundary to the lower boundary. 
This total process is repeated for each time step from m = 1 to the maximum time 
step of interest. 

Only slight modification of this method is needed to solve for the six, coupled 


(fi's of equation (A3-17) 
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The technique embodied in equations (A3-7) (A3-17) and equations (2-25) through 
(2-33) was programmed for the IBM 360-75 for both sets of finite difference equations 
(A3- 7) and (A3-17). The program to solve for 0 was called BOLTSOL and the 
program to solve for the <£. *s was called BOLTDEV. Both programs are listed in 
Appendix IV. Details of the numerical solutions are discussed in Chapter III. 




Figure 2-1. Rates of production and loss of vibrational energy and diffusion 


CHAPTER HI. APPLICATION TO SAR-arc 


The finite difference equations (A3-7) and (A3-17) were solved in single pre- 
cision over the altitude interval of 120 km to 920 km for conditions corresponding 
to those of a 17.7 KR SAR-arc, as shown in Figure 1-2 and 1-3. The upper altitude 
limit was lowered to 840 km. with negligible effect (less than 5° in the vibrational 
temperature and less than 2% in the 0i's) on the solutions. The altitude step was 
taken as 40 km. and variation of the time step by a factor of 10, from 100 sec. to 10 
sec. had no significant effect on the results. The solutions were generated over a 
time interval of 100000 seconds and approached (or reached) a steady state during 
that time. 

Figure 3-1 shows the number of density of vibrational quanta versus altitude 
for a time into the solution of 95200 seconds. The number density of quanta was 
calculated from the solutions of equation (A3-7) corresponding to equation (2-23) 
(that is ignoring the e./n term in the source term). Also shown in Figure 3-1 is the 
quantity i e. calculated from the solutions of equations (A3-17). Theoretically 
equations (A2-4) and(A2-5) should be satisfied, however, because of the finite dif- 
ference and other approximations, i e. will not be zero. A measure of its sig- 
nificance is a comparison with the number of quanta calculated from the solutions 
of (A3-7). From Figure 3-1 it is seen that . i e. represents a negligible (<1%) 
portion of the energy in the Boltz mann solution for altitudes below 300 km. Between 
300 and 600 km altitude, . i e. increases relative to the quantal number density, 
so that at altitudes of 600 km. and above . i e. represents only about 3% of the 
energy of the Boltzmann solution. These conclusions are valid at all other time 
steps. 
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Figure 3-2 is a plot of | .Z e. /e. | versus altitude for each of the six energy 
levels corresponding to a time into the solution of 94700 seconds. This figure is 
also representative of the results at other times. It is seen in Figure 3-2 that be- 
low 400 km altitude . e. is comparable to or greater than the number densities 
of each of the deviations for all 6 energy levels. Above 400 km .2 e. is only 
several tens of percent of the deviations for levels 2 through 5 and a few percent 
of the deviations for levels 0 and 1. It should be pointed out that at 160 km. the 
ratio of the magnitudes of the deviations of level zero to level five (that is € 0 /e s ) 
is 10 7 initially, and for the time corresponding to Figure 3-2, is near 10 6 . This 
explains why the ratios at 160 km. increase with level number. At 840 km. altitude, 
for the time corresponding to Figure 3-2, the ratio of e Q /e g is near 10 2 . It is sig- 
nificant that the ratio | • I e. /e . I has its largest values at the lowest altitudes 

i = 0 1 1 

where the deviations e. are numerically the largest. At low altitudes I e. is 

1 i = 0 1 

about 10 4 of the molecular nitrogen number density and at high altitudes it is about 
10 “ 3 of the molecular nitrogen number density. Figure 3-2 indicates that caution 
should be exercised in using the calculated deviations, e . , only for altitudes below 
370 km. However, as later sections show, the influence of the deviations below 
370 km. altitude is not large, but the e. do become physically important at higher 
altitudes where the precision of determination is considerably improved. 

3.1 Solution for the Boltzmann part. 

The vibrational temperatures calculated from the solution of the \p equation 
are shown versus altitude and time in Figures 3-3 and 3-4, respectively. The 
solutions are for the set 2 (initial and lower boundary vibrational temperature 
of 355°K) boundary conditions. 
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Figure 3-3 shows the vibrational temperature profiles, calculated at five time 
steps, versus altitude. It is seen from this figure that at altitudes above 600 km 
the vibrational temperature is nearly isothermal at each time step, a result of the 
rapidity of the diffusion process at these altitudes. The downward movement of the 
energy is apparent in the figure. Vertical vibrational temperature gradients of as 
much as 25°K/km. are observed in the steady state. 

Figure 3-4 is a plot of vibrational temperature versus time at ten altitudes of 
interest. It is seen from the figure that a steady state is reached in about one day. 
However, in 10 hours, or less depending on altitude, the vibrational temperature has 
achieved 90% or more of the steady state value. Notice that the vibrational temp- 
erature at the higher altitudes has a very steep gradient with time, initially. This 
rapid heating rate is important to the solutions of the deviation equations discussed 
in the next section. It should be noticed that the temperature at 240 km altitude has 
reached 1350°K, a value large enough to significantly effect the reaction constant of 
reaction (1-1), as discussed in section 3.3. 

3.2 Solutions for the deviation parts. 

It was found that convergence of the deviation solutions could not be achieved, 
even with 1 second time steps, if the solution was started at the "turn on" of the 
SAR-arc. It was determined that the difficulty was caused by terms involving 
derivatives of 0 and was apparently due to the initially very large changes in \p 
primarily, those with time. Convergence of the solutions of the ( P^ equations was 
obtained if the solutions were started at 5200 seconds "into" the arc, and this was 
the approach adopted. The physical conditions corresponding to the ini tial stages 
of a SAR-arc are unknown. However, it is expected that the characteristic times 
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of the important processes are measured in terms of one or more hours. Thus, 
the imposition of SAR-arc conditions in a step-wise manner is not realistic; so 
that little physical content has been lost by the adopted calculational approach. 

Figures 3-5 through 3-10 show the <£.'8 versus altitude for three time steps. 
These figures show that above 500 km the ratio of the deviation from the Boltzmann 
solution to the Boltzmann population, at the energy level of interest, is nearly 
constant with altitude. This result is again due to the rapidity of the diffusion 
process at these altitudes. It is also apparent from these figures that at altitudes 
of 300 km and above the 0 . 's rapidly reach a maximum and then decrease to their 
"steady state" values. In these figures, where 4> i changes sign with altitude the 
absolute value has been plotted. It is seen that the "steady state" deviation popu- 
lations above 300 km altitude, are less than the Boltzmann population of the cor- 
responding energy levels 0, 1, 2 and 3 and are greater than the Boltzmann popula- 
tion for levels 4 and 5. This will be important in the calculation of an effective 
reaction coefficient for reaction (1-1) (see section 3.3). 

Figures 3-11 through 3-16 show the 0. 's versus times for several altitudes. 
From these figures it is seen that above 300 km altitude a nearly constant, sig- 
nificant, deviation is maintained from a Boltzmann distributed population for the 
upper energy levels. 

3.3 Reaction rate constant for reaction (1-1) and ionospheric loss time constant. 

The effective reaction rate constant for the ion-atom interchange reaction (1-1) 


was calculated from 
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i = 0 



K + g i) 

n 


(3-1) 


where k e is the effective rate constant and k(i) is the rate constant for level i, taken 
from reference (13) and listed in Table 3. The reaction rate for the Boltzmann part 
of the solution, kg, was calculated from equation (3-1) with the e. terms set to zero. 
The computer program for these calculations, REACTRAT, is listed in Appendix HI. 

A comparison of k e and k B is a measure of the importance of the deviation from a 
Boltzmann distribution to the F 2 region ionospheric ionization loss. 

Figure 3-17 shows the ratio k g /k B versus time for several altitudes of inter- 
est. It is apparent in the figure, that for altitudes above 320 km the ratio is sig- 
nificantly larger than 1 throughout the time interval of the solutions, and particularly 
during the first 10 hours, the approximate duration of SAR-arcs. 

Figure 3-18 shows the ratio k e /k B versus altitude for five time steps. It is 
seen in the figure that for a given time step the ratio is relatively constant with 
altitude for altitudes above 500 km. There is a large variation of the ratio with 
altitude between 280 km and 500 km, the region of the atmosphere where diffusion 
becomes the dominant process. This figure shows that the deviations from a Boltz- 
mann solution are noticeable for altitudes as low as 300 km. 

Figures 3-19 and 3-20 show the effective rate constant, calculated from equa- 
tion 3-1 using the ^ and 0. equation solutions, versus time and altitude, respectively. 
These figures show that the rate constant can be factors of 2 to 8 times larger than 


k e/ k 
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the value for vibrational temperatures near 300°K (shown as a triangle in the two 
figures). 

The k e values, shown in Figures 3-19 and 3-20, have been used to calculate 
the loss rate for the ionospheric electrons and the ionospheric 'loss time constant" 
in the following way. Assuming reaction (1-2) is infinitely fast, then the electron 
loss rate will be the same as the 0 + ion loss rate. Further assuming the electron 
density and the 0 + ion density are the same, we have 



= - k e n n e 


(3-2) 


where ne is the electron density. The loss rate for the electrons is defined from 
(3-2) to be 


— = k n (3-3) 

r e 

where t is the ionospheric loss time constant. 

Figure 3-21 and 3-22 show the loss time constant calculated from (3-3) versus 
altitude and time respectively. In Figure 3-21 the loss rate for four time steps and 
the "normal" (vibrational temperature of 300°K) ionospheric time constant are 
shown. Noticeable decreases in the ionospheric loss time constant from the "nor- 
mal" value is observed at altitudes as low as 280 km. 

Figure 3-22 shows the loss time constant versus time for several altitudes 
where the time constant is one day or less. It is seen from the figure that in the 
first 10 hours, the time constant at 360 km. altitude, where the "normal" value is 
approximately 90 minutes, decreases to a value of 30 minutes. This figure shows 
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that the ionospheric loss time constant is decreased by factors larger than three 
in altitude regions between 360 km and 440 km where the ionospheric loss time 
constant is measured in hours. 

3.4 The effects of: neglect of e, /n; boundary conditions; and atomic oxygen loss 
process. 

To determine the effect of the approximation of neglecting the e. /n term in the 
right hand side of the 0 equation, (A2-30), the solutions of the <£. equations were 
used to calculate the e /n terms. Equation (A2-30) was then solved using those 
e. /n terms. The second approximation to </> was then used to recalculate the 4> l t s. 

The second approximation to ^ yielded a vibrational temperature which was 
different from that obtained from the first approximation by 67°K, or less, at higher 
altitudes. At lower altitudes the difference was only a few degrees. The I i e 

i 3 0 * 

terms from the new solutions were also reduced by as much as a factor of 3. Thus, 
as expected on the basis of the . 2 q i e. results from the first solution, the second 
approximation to <// produces only a few percent correction to the results of the 
first approximation. The e. 's calculated from the second approximation differed 
from those of the first approximation by 30%, or less, with most of the difference 
occurring at altitudes below 360 km. The . e. term changed by a few tens of 
percent only. It is concluded, that the second approximation produces a negligible 
correction to the 4 1 equation. Because most of the influence of the e 's is at alti- 
tudes above 300 km. (see section 3.3), the second approximation to the e . ’s is 
different from the first by only several percent in its impact on the ionospheric 
electron loss rates, and is therefore also negligible. 
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The and 4> i equations were resolved using set 1 boundary conditions. It 
was found that the vibrational temperature calculated from the 0 equation solution 
differed from that obtained using set 2 boundary conditions by only a few degrees . 
The €. results differed by less than 3%. These results indicate that the solutions 
are not sensitive to the details of the loss mec hani sms at the lower boundary. 

An interesting result was obtained by comparing the behavior of the solutions 
with and without inclusion of the atomic oxygen production and loss processes for 
vibrational excitation, for set 1 boundary conditions (that is T v (120, t) = T v 
(z, 0) = 0°K). The solution including the oxygen process yielded vibrational tem- 
peratures that were larger at all altitudes than those obtained from solutions not 
including the atomic oxygen processes. The largest differences between the two 
solutions initially occurred at altitudes below 320 km and was greater than 150°K. 
This location for the differences is expected because this is the altitude region 
where the oxygen effects are relatively the most important. This result can be 
explained by the fact that the boundary condition of 0°K vibrational temperature 
causes the initial vibrational temperatures to be less than the gas temperature, 
thus the atomic oxygen is a source of vibrational excitation through translational- 
vibrational energy exchange. Eventually the electronic source of vibrational 
energy dominates the solution. In the "steady state" the vibrational temperature 
corresponding to the solution with no atomic oxygen processes has a vibrational 
temperature only 8°K greater at 160 km altitude and 6°K greater at 920 km alti- 
tude than that corresponding to the solution containing the atomic oxygen processes. 

A comparison of the <f> equation solutions for these two cases is complex. At 
altitudes below 280 km, the e. initially changed by factors as large as 30 and 40 
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with the largest values of e. corresponding to the solutions containing the oxygen 
processes. Above 320 km altitude the two solutions for the e. »s differ by ap- 
proximately 25% or less. The e. x *s of the upper two energy levels show the 
smallest differences between the two solutions. The two sets of solutions for the 
e. 's approach the same values in the steady state. The influence of the differ- 
ences in the two sets of solutions on the ionospheric loss rates is small because 
the large differences occur below 320 km. altitude where the deviations from the 
Boltzmann part of the solutions play a small part in the determination of the 
reaction rate constant. 



CHAPTER IV. SUMMARY AND CONCLUSIONS 


The objective of this work was to determine if the observed decrease in 
ionospheric electron density occurring in SAR-arcs could be explained by an in- 
crease in the reaction constant for the rate determining ion-atom interc han ge re- 
action (1-1). The increase in the reaction constant is due to vibrational excitation 
of molecular nitrogen by subexcitation, inelastic collisions between ambient, 
thermal, electrons and nitrogen molecules. The finding of this work is that the 
vibrational excitation of molecular nitrogen is indeed important to the ionosphere 
under SAR-arc conditions and can cause the loss time constant for the F 2 region 
to change to lower values by factors of 1.6 to 5 with an attendant decrease in the 
electron density. 

Previous work, (20), (22) on the related problem of the diur na l variation of 
the molecular nitrogen vibrational temperature assumed the vibrational distribu- 
tion to be Boltzmann and considered diffusion processes. Earlier work (23) in the 
vertical distribution of vibrationally excited molecular nitrogen in aurora ignored 
diffusion and solved the continuity equations for the excited energy levels of 
molecular nitrogen. This work solves the continuity equations for the first six 
energy levels of molecular nitrogen including the diffusion process. Thus, it 
combines into one solution, the previous' approaches. 

The continuity equations to be solved are six, coupled, mildly non-linear, 

partial differential equations contai ning the dominant source and loss terms for 

vibrational excitation of molecular nitrogen. The solutions were divided into an 

energetically equivalent part that is Boltzmann distributed and a part representing 

deviations from the Boltzmann solution. The resulting equations were solved 
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numerically, using finite difference equations, to provide the populations of the 
energy levels as functions of time and altitude. 

It was found, under conditions where the ionospheric electron temperature 
reaches values near 5000 °K at 500 km altitude, that the vibrational temperature 
of molecular nitrogen reached a steady state value near 2300°K. at that altitude. 
The deviations from the Boltzmann distribution depends strongly on altitude for 
altitudes below 400 km. This is the region of the atmosphere where the diffusion 
process becomes dominant with increasing altitude over the vibrational-vibrational 
energy exchange process. It was found that for the highest three energy levels, 
above 300 km. altitude, the deviations from a Boltzmann distribution increased 
the populations of these levels by factors of 2 to 10. These results were insensi- 
tive to plausible variations of the boundary conditions. 

The deviations from a Boltzmann distribution caused the effective reaction 
rate for reaction (1-1) to be increased by as much as a factor of 1.7 over the value 
obtained from the energetically equivalent Boltzmann solution. This indicates the 
possible importance of these deviations in problems involving chemical reactions 
at altitudes as low as 500 km. The effective reaction rate, calculated using the 
total population densities, was as much as a factor of 8 greater than the reaction 
rate for a vibrational temperature of 300°K. Significant enhancements of the re- 
action rate were found to occur at altitudes as low as 300 km. 

Because of these increased reaction rates, the times required to significantly 
decrease the electron density in the ionosphere were reduced from their normal 
value by factors of 1.6 to 5, depending on altitude. Over most of the altitude range 
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where the loss time is short compared to SAR-arc durations, the loss time was 
reduced by factors exceeding 1.8. 

These results lead to the conclusions that the vibrational excitation of molec- 
ular nitrogen, under conditions corresponding to SAR-arcs, causes a large increase 
in the loss rate of ionospheric electrons in the F 2 region. The enhanced loss rate 
is of such magnitude, .that it alone will account for the observed electron density 
depressions in SAR-arcs. 



REFERENCES 


1. Barbier, D., "L'Activite aurorale aux basses latitudes," Ann. Geophys., 14, 
334, 1958. 

O 

2. Roach, F. E. and Roach, J. R., "Stable 6300 A auroral arcs in mid-latitudes," 
Planet. Space Sci., 11, 523, 1963. 

3. Nagy, A. F., Roble, R. G., Hays, P. B., "Stable mid-latitude red arcs: ob- 
servations and theory," Space Science Rev., 11, 709, 1970. 

4. Cole, K. D., "Stable auroral red arcs; sinks for energy of D st main phase," 
J. Geophys. Res., 70, 1689, 1965. 

5. Roble, R. G., Hays, P. B., and Nagy, A. F., "Comparison of calculated and 
observed features of a stable midlatitude red arc," J. Geophys. Res., 75, 
4261, 1970. 

6. Roble, R. G„ Norton, R. B., Findlay, J. A., and Marovich, E., "Calculated 
and observed features of stable auroral red arcs during three geomagnetic 
storms," J. Geophys. Res., 76, 7648, 1971. 

7. Chandra, S., Maier, E. J., Troy, Jr., B. E. and Narasinga Rao, B. C., "Sub- 
auroral red arcs and associated ionospheric phenomena," J. Geophys. Res., 

i 

76, 920, 1971. 

8. Nagy, A. F., Hanson, W. B., Hoch, R. J., Aggson, T. L., "Satellite and 
ground-based observations of a red arc," J. Geophys. Res., 77, 3613, 1972. 

9. Norton, R. B. and Findlay, J. A., "Electron Density and temperature in the 
vicinity of the 29 September 1967 middle latitude red arc," Planet. Space 
Sci., 17, 1867, 1969. 


64 



65 


10. Roble, R. G., "A theoretical and experimental study of the stable mid- 
latitude red arc (SAR-ARC)," Ph.D. Thesis, University of Michigan, 1969. 

11. Rishbeth, H. and Garriott, O. K., Introduction to Ionospheric Physics , 106, 
Academic Press, New York, 1969. 

12. Schmeltekopf, A. L., Fehsenfeld, F. C., Gilman, G. I. and Ferguson, E. E., 
"Reaction of atomic oxygen ions with vibrationally excited nitrogen mole- 
cules," Planet, Space Sci., 15, 401, 1967. 

13. Schmeltekopf, A. L., Ferguson, E. E. and Fehsenfeld, F. C., "Afterglow 
studies of the reactions He + , He (2 3 S), and 0 + with vibrationally excited N 2 ," 

J. Chem. Phys. 48, 2966, 1968. 

14. Jacchia, L. G., "The temperature above the thermopause," Space Res. , 5, 

1152, 1965. 

15. Dalgamo, A., "Vibrationally excited molecules in atmospheric reactions," 
Planet. Space Sci., 10, 19, 1963. 

16. Dalgarno, A., Me Elroy, M. B., and Moffett, R. J., "Electron temperatures 
in the ionosphere," Planet. Space Sci., 11, 463, 1963. 

17. Hunten, D. M., "Some photometric observations of auroral spectra," J. 

Atmos. Terr. Phys., 7, 141, 1955. 

18. Thomas, L. and Norton, R. B., "Possible importance of internal excitation 
in ion-molecule reactions in the F region," J. Geophys. Res., 71, 227, 1966. 

19. Walker, J.C.G., "Electron and nitrogen vibrational temperature in the E- 
region of the ionosphere," Planet. Space Sci., 16, 321, 1968. 

20. Walker, J.C.G., Stolarsky, R. S. and Nagy, A. F., "The vibrational temperature 
of molecular nitrogen in the thermosphere," Ann. Geophys., 25, 831, 1969. 



66 


21. McNeal, R. J., Whitson, Jr., M. E. and Cook, G. R., "Quenching of vibration- 
ally excited N 2 by atomic oxygen," The Aerospace Corporation SPL 3303, 

L. A. Cal., July 1972. 

22. Breig, E. L., Brennan, M. E. and McNeal, R. J., "Effect of atomic oxygen 
upon the N 2 vibrational temperature in the lower thermosphere," the Aero- 
space Corporation SPL Report 3301, L. A. Cal, July 1972. 

23. Schunk, R. W. and Hays, P. B., "Theoretical N 2 Vibrational Distribution in 
an Aurora," Planet. Space Sci., 19, 1457, 1971. 

24. Walker, J.C.G., "Oxygen and Nitrogen Vibration in the Thermosphere," 
paper presented at Summer Advanced Study Institute, Physics and Chemistry 
of Upper Atmosphere, Univ. of Orleans, Orleans, France, 9 August 1972. 

25. Kummler, R. H. and Bortner, M. H., "Vibrational temperatures in the E and 
F regions," Space Research XII, 711, Akademie-Verlag, Berlin, West 
Germany, 1972. 

26. Bauer, E., Kummler, R. and Bortner, M. H., "Internal energy balance and 
energy transfer in the lower thermosphere," App. Optics, 10, 1861, 1971. 

27. Herzfeld, K. F. and Litovitz, T. A., Absorption and dispersion of ultrasonic 
waves, Academic Press, New York, 1959. 

28. Chen, J.C.Y., "Theory of subexcitation electron scattering by molecules. I. 
formalism and the compound negative ion states," J. Chem. Phys., 40, 

3507, 1964. 

29. Chen, J.C.Y., "Theory of subexcitation electron scattering by molecules. II. 
excitation and de-excitation of molecular vibration," J. Chem. Phys., 40, 


3513, 1964. 



67 


30. Chen, J.C.Y., "Erratum: theory of subexcitation electron scattering by mole- 
cules. n. excitation and de-excitation of molecular vibration," J. Chem. Phys., 
41, 1964. 

31. Treanor, C. E., Rich, J. W. and Rehm, R. G., "Vibrational relaxation of an- 
harmonic oscillators with exchange-dominated collisions," J. Chem. Phys., 

48, 1798, 1968. 

32. Fisher, E. R. and Kummler, R. H., "Relaxation by vibration- vibration exchange 
processes, part 1, pure gas case," J. Chem. Phys., 49, 1075, 1968. 

33. Sharp, F. E. and Rapp, D., "Evaluation of approximations used in the calcu- 
lation of excitation by collision. I. vibrational excitation of molecules," J. 

Chem. Phys., 43, 1233, 1965. 

34. Hirschfelder, J. O., Curtiss, C. F. Bird, R. B., Molecular theory of gases and 
liquids, John Wiley and Sons, New York, 1954. 

35. Chapman, S. and Cowling, T. G., "The mathematical theory of non-uniform gases, 
Cambridge Univ. Press, London, 2nd edition, 1960. 

36. Richtmyer, R. D., Difference methods for initial-value problems, Interscience 
publishers, Inc., New York, 1957. 

37. Feshback, H., "Unified theory of nuclear reactions," Ann. Phys., 5, 357, 1958. 

38. Davydov, A. S., Quantum Mechanics, trans. by D. ter Haar, Pergamon press, 


New York, 1965. 



68 


TABLE 1 

Cross Sections used from Reference (30) in Equation (Al-19) 


Cross-Section 

Comment 

CT 10 

Numerical data available 

a 20 

Numerical data available 

°03 

Numerical data available 

CT 04 

Numerical data available 

°- 5 o 

All that is available 

a !2 

All that is available 

CT i3 

All that is available 

CT 14 

All that is available 

CT S1 

All that is available 

CT 23 

All that is available 

CT 24 

All that is available 

a 52 

All that is available 

Cr 34 

All that is available 

CT S3 

All that is available 

CT S4 

All that is available 
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TABLE 2 

Parameters for Reaction Rate Constants for Production and 
Loss of Vibrational Excitation by Electron Collisions 


I.D. 

PI 

P2 

P3 

P4 

Max Error % 

10 

-1.837483E01 

-2.218530E04 

1.684328E07 

-4.558864E09 

-19@2500° 

20 

-1.887149E01 

-1.851033E04 

-1.31 7388 E06 

2.310335E09 

-37@1500° 

50 

-2.004471 E01 

-1.797215E04 

-6.671288E06 

2.076589E09 

-14@2000° 

51 

-2.025835E01 

-1.628830E04 

-2.509007E06 

7.790528E08 

- 5@1500° 

12 

-1.874858E01 

-1.222789E04 

1.102594E07 

-2.9931 75E09 

14@5000° 

52 

-1.974898E01 

-1.314905E04 

-2.739562E06 

9.167196 E08 

- 9@1500° 

03 

-1 .950456 E01 

-5.705000E03 

-6 .362043 E06 

1.91 7276 E09 

-13@1500° 

13 

-1.972427E01 

-4.652598E03 

-6 .585213 E06 

1.93141 9 E09 

-12@1500° 

23 

-1.939743E01 

-8.030383 E03 

5.607227E06 

-1.257535E09 

-16@2000° 

53 

-1.980759E01 

-8.800227E03 

-2.526216 E06 

9.740170E08 

-12@1500° 

04 

-1.964117E01 

-3.913067E03 

-5.735326 E06 

1.842685 E09 

-13@1500° 

14 

-2.019481E01 

-3.588698E03 

-4.862049 E06 

1.688892E09 

- 9@1500° 

24 

-1.999716E01 

-4.4211 76 E03 

-3.194914E06 

1.243446 E09 

-10@1500° 

35 

-1. 943835 E01 

-4.869426 E03 

2.078704E06 

-1.935087E08 

-15@1500° 

54 

-1.983298E01 

-4.700520 E03 

-4.21 7622 E05 

4.723745E08 

-14@1500° 




TABLE 3 


Reaction Rates for Vibrational Levels 


Level 

k y [cm. 3 sec. -1 ] 

0 

1.3 

x 10“ 12 

1 

1.3 

x 10“ 12 

2 

5.0 

x 10" u 

3 

1.20 

x 10‘ 10 

4 

2.6 

x 10 _1 ° 

5 

3.6 

x 10" 10 

6 

4.4 

x 10 -10 

7 

4.6 

x 10" 10 

8 

4.6 

x 10 -10 

9 

4.6 

x 10 -10 

10 

4.6 

x 10 -10 

11 

4.6 

x 10“ 10 






APPENDIX I 


ELECTRONIC PRODUCTION AND LOSS OF VIBRATIONAL EXCITATION 

Production and loss rates for the first six vibrational levels of molec ular ni- 
trogen were calculated using theoretically calculated (28), (29), (30) cross-sections 
for inelastic, sub-excitation, electron scattering for these levels. The cross- 
sections were averaged over the Maxwell- Boltzmann energy distribution repre- 
senting the electron gas to produce the production and loss rates. This section 
contains an amplified outline of the content of the theoretical cross-section calcu- 
lations in Chen's, (28), (29) papers as well as discussion of the calculation of re- 
action rates. 

Chen (28), (29) applies the unified theory of Feshbach (37), for nuclear reac- 
tions involving the formation of a compound state, to inelastic scattering of elec- 
trons by molecules. Herafter Chen (28) will be called paper I and equations con- 
tained therein referenced as I X.XX. Chen (29) will be referenced as paper II and 
equations contained therein will be referenced as n X.XX. 

A1.1 Paper I: 

Paper I starts with the consideration of an N-electron system consisting of an 
incident electron and a diatomic target molecule with N-l electrons. The system 
is separated into relative and center of mass coordinates and the center of mass 
part is removed from the Schrodinger equation. The relative coordinates for the 
electrons, X, , ? 2 . . . , ? N-1 , are referenced to the center of mass of the 
nuclei. X refers to the incident electron. R is the vector separation of the two 
heavy nuclei a and b with masses M a and M fa respectively. The Schrodinger 
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equation is (in atomic units) 


[K 0 (X) + V 0 (X, 7, R) + H ab (r, R)] ft = Eft 


( 12 . 1 ) 


where 



( 1 - 2 . 2 ) 


and r is the collective notation of r x , r 2 , . . . , r N _ x for the molecular electrons; K 0 
is the kinetic energy operator of the incident electron; V 0 is the Coulombic poten- 
tial energy of the incident electron; H ab is the Hamiltonian for the target molecule; 
ft is the total wavefunction for the system having the total energy E; X a and X fa are 
the position vectors of the incident electron with respect to the two nuclei having 
and Z b as their charges. 

It is assumed that the unperturbed molecule is described by the set of 
Schrodinger equations 


where 


and 


H el (7, R) $(r, R) = e n (R) ® n (r, 5) (1-2.3) 

[- (2m) _1 % + fi n (R)] X nr (R) = e ny Xny (R) (1-2.4) 


«,(?, R) = 


1 + M a + \ 

2 (M a + M,) 


N-l 

22 w +«.,<?. r) 

i=i 


(1-2.5) 


u2 

M = ; V: 

M a + ^ 1 

is the Laplacian operator on coordinate i and U el 


is the potential energy of the 


molecular electrons 
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H ab (7, R) = - (2/x)' 1 V* h- H el (r, R) (1-2.6) 

The cross terms for electrons and electron-nuclei contained in the Born-Oppen- 
heimer separation are neglected because they do not enter, in the zeroth order ap- 
proximation, in the expression for the excitation of nuclear motion, (34), pages 
925-928. The electronic wavefunctions, $ n , and nuclear wavefunctions y ny hav- 
ing eigenvalues e n (R) and e ny , respectively, form a complete, orthogonal set in 
their corresponding spaces. 

The total wavefunction is expanded in terms of <$ n and y to give 

fi(X, 7, R) = 22 R) y nr (R) 0 n7 (X) (1-2.7) 

n,r 

where exchange effects are neglected. The expansion coefficients </* (x) are as- 

sumed not to depend on the inter-nuclear separation as a parameter. It is as- 
sumed that the inter-nuclear separation is fixed at its equilibrium value R( 0 ) for 
0 ny (X). The effect of this assumption is estimated to be of order /T 1 . 

Substituting equation (1-2.7) into the Schrodinger equation (1-2.1), multiplying 
on the left by $*,, \* n , y , , integrating over r and R and using (1-2.3) through 
(1-2.6) we have the coupled sets of equations 

<*.<*) + - El ^ *.'■>" ( *> 

n' ,r' 

with 

V ”' (X) = (X nr IV„„.(X. 5 ) - (2m)-‘ C„„.(R, v B )| Xn , y ,) 
v„„,(x, R) = <®|v 0 (X, R, ?)!*„,) 


( 1 - 2 . 8 ) 


(1-2.9) 


( 1 - 2 . 10 ) 
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c„„. = 2 -7, + <*„|vllv d- 2 - 11 ) 

Because we are interested in subexcitation electron impact, the V' 0 -y’ s are of par- 
ticular interest. Thus writing (1-2.8) as 


M) + - a -- - £ v 0 «" v * . £ 0 nv( it) 

n' ,y' 
n'/o 


(Al-1) 


y‘ 


and 


1*0 <*> + e ny - E] 



for n/0. 



(Al-2) 


Rewriting (1) and (2) we have 


[K 0 (X) + e 0y ] ^(X) + ^ Vft" 0^, (X) - E^(X) = _ ^ ® 0 n , yf (X) (Al-3) 

y' n'y' 

n'j'O 

+ e »y ] + £ V ”' (X) - E *nrW * -£ V"' ^ (X) (Al-4) 

n'.y' y' 

n'/O 


Defining 


+ e 


nO 


T = 

n 


K 0 + e nl ® 


^0 + e n2 


(1-2.14) 


0 
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( 1 - 2 . 12 ) 


( 1 - 2 . 15 ) 


( 1 - 2 . 13 ) 


( 1 - 2 . 16 ) 


( 1 - 2 . 17 ) 
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(3) and (4) become 


T 0 ^0 + U 00 ^0 - E h = - 22 U 0n ' € = - U 0 * 


and 


n* 70 


or 


T \b + V U , dj , - E0 =-U„0 A 

n r n / nn' r n' r n nO r 0 


n'/O 


where 


(H - E) V = - u; 0 O 



T 1 + U H 

U 12 

U 13 



u 21 

^2 + ^22 

U 23 

* 

H = 

U 31 

U 32 

T 3 + U 33 

* 


and 


E f = E8 , 

nn' nn' 


Following Fishbach’s formalism, (37), the formal solution of (1-2.20) is 


V = (E - H + i v y l u; 0 O 

where 77 = 77 8. . and 77 o + to ensure outgoing waves in the exit channels, 
gives 0 in terms of $ Q . Using (1-2.21) in (1-2.19), we obtain 


(1-2.19) 


( 1 - 2 . 20 ) 

(1-2.18) 


( 1 - 2 . 21 ) 

( 1 - 2 . 21 ) 


(T 0 + A - E) 0 O = 0 


( 1 - 2 . 22 ) 
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where 

A = U 00 + u 0 (E - H + iri y l u; (1-2.23) 

and is the effective scattering potential for subexcitation scattering of electrons by 
molecules. 

The formal solution (1-2.21) contains the Green's function operator lim 

n-0 + 

(Ej - H + i^ 7 )“ 1 which implies the need for a set of eigenfunctions of the effective 
Hamiltonian H. Thus we write 


«’*'»= (1-3.2) 

where a set of eigenfunctions {<P m } approximates the set of compound negative ion 
eigenfunctions with energy spectrum {W } . 

m 

The physical significance of these states is illustrated in the weak coupling 
limit where H is diagonal. In this case (1-3.2) gives the equations 

( T n + U n„> ^ = »<"> 'A!" 0 " / 0 (1-3.3) 

which is expanded further to become 

»o * V" + - W<5>] -- - £ (1-3.4) 

y'fy 

The compound negative ion states <//<>> describe, in this instance, a set of states 
of the incident electron in resonance with the field of the excited target molecule. 
Expanding the Green's function operator contained in (1-2.21) in terms of the 


eigenfunctions l t l m we have 
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g (+) (e> = y 


I'M 

m n 

E - W 


•J 


dp dW 


I»(P. W )><¥(p, W)| 

E - W + ip 


(1-4.1) 


where the sum over m is for the discrete states and the integral over p contains the 
continuous eigenspectrum for energy E. Expression (1-4.1) implies a resonance part 
of the operator and a non-resonance part. Therefore we write 


G< + >(E) = Q(E) + G< + >(E) 


(1-4.2) 


with 




(1-4.3) 


Q(E) is the resonance part of the operator and the prime means that only those 
negative ion states having energies comparable to that of the system are included 
in the summation. The effective potential (1-2.23) is now written, using (1-4.2) in 
(1-2.23) as 


where 


A = u* + u 


(1-4.4) 


u' = U 0 Q(E) U 


o 


(1-4.5) 


u=U oo + u o G P( £ ) % 


Let 77 be the solution for an outgoing wave of the equation 


(1-4.6) 


(T 0 - E) tt(-) r 0 
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Using (1-4.4) in (1.22) we have 

(T 0 + U-E)0 o = -u' 0 O (Al-5) 

where the resonance part of the effective scattering potential is treated as the 
scattering potential. 

From scattering theory, (38), Chapter 12, the solution for the outgoing wave is 
to = ^ +) + (E - T 0 — u + i-77)- 1 u' 0 O (1-4.12) 

where 

(T 0 + u - E) 0 o < + > = 0 (1-4.10) 

and + ) is an incoming wave and contains the effects of the direct scattering 
potential u. 

The transition matrix describing transitions from the incident channel y' to 
channel y for n = 0 is 


3 0 (% V) = |u +u'| 0 O ( y')> (1-4.13) 

Consider the outgoing wave solutions to 

(T 0 +u -E) 0 o <-> = 0 (Al-6) 

which may be written as 


0 o (-) = 7T 0~ ) + ( E - T 0 - i7 7) -1 u ^o~ ) 


(1-4.14) 
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Multiplying by (E - T 0 - i 17 ) and adding and subtracting U 7 r(“> to the right hand 
side we have 


(E - T o - 177 ) «//<-> = (E - T 0 - i 77) ^ 0 (_) + + U77 0 <~) -U 77 0 C> (Al-7) 

>/'o ( " ) = 77 o" ) + ( E - T o - u - i 7 ?)" 1 u 77 o" ) (1-4.15) 

- "P = 'I’P - (E - T 0 - u - iT,) u 77 0 (-> (Al- 8 ) 

Substituting (A-l- 8 ) in (1-4.13) we have 


3 0 (% y') = < 7 T o~ ) l u l </'o> - <( E - T o - u - i7 ?) -1 u 77 o (_) l u ' I 'V 

+ < l /' 0 ( " ) |u'| 0 O > 


(1-4.16) 


= <77- 0 (") |u [ 0 O - (E - T 0 - u + iT?) -1 u' I p 0 > 

+ < |u' | sA 0 > (Al-9) 

= < 77 0 (-) W 0^ +) > + <'A 0 (-) |u'| <// Q > (1-4.17) 

using (1-4.12). Consider |u* | t// 0 > and repeatedly substitute (1-4.12) for </> 0 
to obtain 

< l /' 0 ( " ) |u'| 0 O >= <0 O (_) lu'l 0 o (+) > +< 0 o (-) |u'(E - T 0 - U + u' |t p 0 ) 

= l u> I ^o ( + ) > + l u * (E - T 0 - u + irj)~ l u ' |i/£ +) > 

+ |u'(E-T 0 -u + It ,)" 1 u(E-T 0 -u + i7 ? )- 1 u' |</'o> 

+ • • * • 

= <0 O <’> |u' +u' (E - T 0 - u + i 77 ) _i u' + u'(E-T 0 -u + i7 7 )- 1 u' 

‘ (E - T 0 - u + i-ry) -1 u' + |0 o ( + ) } 


(Al-10) 
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For an isolated resonance, the summation over m in (1-4.3) reduces to one term 
and (1-4.5) becomes 


Thus (11) becomes 


u = 


Kx'f. K 

E - W 


(Al-11) 


Let 


<>/' 0 ( " ) |u' I 0 O > 


<^P lUol’.XV. K\ H + ) > 

E - W 

m 


f <f. IU;(E-T 0 -U t i,)-' U 0 »„) 

I' ^ 



K = <*. lu; (E - T. - u + i v T' U 0 |»> 


Therefore 


<'A 0 ( ' ) l u ' I 0 O > = 


<C > l u 0 l *«><*„■ l u ol ^ 0 ( + ) > 
E - W 


R 


E - W E - W 


f U ol ^ +) > 


E - W - R 

m m 


(1-4.23) 


(1-4.24) 


A1 .2 Paper n. 

In his application to molecular nitrogen, (29), Chen makes the following 
assumptions: 

(1) an isolated resonance occurs. 

(2) the direct transition amplitude is negligible. 
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(3) a single excited electronic state, n, of the target molecule is responsible 
for the resonance. 

(4) the nuclear states of the target molecule are not strongly coupled and the 
Born-Oppenheimer coupling matrix elements are negligible. 

(5) the total contribution to the transition amplitudes due to resonance effects 
of the set of compound negative-ion states may be accounted for as a linear 
sum of the individual contribution of each compound state. 

From assumptions (1) and (2) equation (1-4.24) gives the transition matrix ele- 
ment. Thus 


o 0 (% Y) = 


<V'o ( ' ) (7) |U 0 I lU'l *o ( + > </)> 


E - W - R 

m m 


(n-2.2) 


From assumption (3) only V on , V no , V nn and V QO in equation (1-2.10) are non- 
zero. Thus H in equation (1-2.18) becomes a single element matrix. Assumption 
(4) gives from equations (1-2.11) and (1-2.9) respectively 


C, = 0 and V”" = 0 - 

nn nn 


Thus U nn in equation (1-2.13) is diagonal. Therefore we have 


giving 



H = W (m > 

n r n n r n 

or 


(T + U ) = W (m) 

v n T nn / ' n 


(Al-12) 



where 
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0 

0 


0 n (m) (M 



0 


Expanding (Al-12) further and rearranging terms 




n X 


X 


'nO' 


"nX 


Now consider 


Similarly 


Jy 

= (u„„ ^ <m, V 

= *<■> 



(u;* 


|0 ) - V y '^ t 

>y on V nk 


(II-2.4) 


(Al-13) 


(Al-14) 


The y - y transition due to the X th vibrational level in the negative ion state is 


3o(% 7') = 


<*&’ Iv^l ^>> (^-x'lv^'l 


n k 1 n0 


or # 


E - W<"> + R<?> 

n A. n A. 
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Using assumption (5), the total transition matrix element is 


o 0 (% /) = ^ 3^ (% /) 


-L 

k 


l V .^ I IV„V' I Q 

E - W<"> + R<!"> 

n A. n A. 


And 


Expanding 


n A. ^ A. ^ L v „V (* - H„ + io)’ 1 V» x !^>) 


lA&^x 'Af I 

'° T ‘'" - 1“^ "E - E + i„ f r 


(E-H. + ir,)' 1 = dE 




where £ is the density of states in the exit channel y , and using 


E - E.. 


177 


E-E r + ii) (E -E y ) 2 +1J ’ <E- Ey > ! + i) ! 


we have 


R (m) . A (n>) _ i J_ r<f > 

nA. nA. 2 n A. 


where 


EE 

A(>") = /!//(■") I V " 1 p f v^' y 2 L 


T? 2 




(Et-2.3) 


(H-2.5) 


(H-2.5) 


(n-2.8) 
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and 

rj™> = 2 Am (R<"> ) 


.= 2 / 0<"> 
n A. 



^IC ) ><C ) l 

(E - E r ) 2 + rf 


^ dI v i *.«*» 


= 2 



IC’>< C’> <*£’ l” 8 ( E - V fy <*Eyl 0^>) 


(Al-15) 

(Al-16) 


Using 


and 


then 


dfl f k y d(k 2 ) 

(2tt) 3 dE y 

dn f ! k y d(k^) 

^32 dE y 


dn. 

iy = — - ky 
(2") 3 


Using (Al-17) to (Al-16) we have 


(Al-17) 


r<?> = — _ 

(2tt) 2 



(H-2.9) 


Using the definitions of and V^y' from equation (1-2.10) and interchanging the 
order of integration in (II-2.3) we obtain 


3 0 <% y 





<^0y[yJ XnAXXHxKr^Xoy 0 
E - E<"> + il/2r<“> 

nA. n A. 


(n-2.il) 



and 
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r »' >= H ^T 2 f W'VJXn^ 8 ^. < n - 2 - 12 ) 

‘~y‘ (27t) j 

where 



(II-2.13a) 


(II-2.13b) 

E «. )= + A<-> = .<»> t AW + ^ 

(H-2.14) 


The differential partial cross section d cr ^, , is expressed as a function of the inci- 
dent angles of the impinging electron by 

dOyy, = (477 2 )" 1 |3 0 (% /)| 2 dfi f 

where 

k 2 = k 2 , + 2(6^, - e 0r ) 

with the initial and final wavenumbers of the electron given by k , and k y 
respectively. 

Because the rotational motion of the molecule is relatively slow, the molecular 
orientation is considered fixed throughout the vibrational excitation or de-excita- 
tion. If the molecule remains in the ground rotational state then the partial cross- 
section can be obtained by averaging the. transition matrix elements over the 
incident angles and integrating the resultant over the sphere of the final an g les. 

The result is 

= (16. 7T 3 )” 1 z/)| 2 dfi 0 dfl f (D-2.17) 



87 


where the y ' s have been replaced by v’s to reflect that the x’s now represent 
vibrational wave functions only. 

The further assumption that the compound negative-ion state does not depend 
on the vibrational state is made. This assumption allows the cross-section to be 
written as 


with 


cr i = 
vv 


yfl B , k 

~v “v V 


1 6 77 3 


L 




(»> = y -1 

n\ 


E - + i — rC 1 ?) 

nA. 2 nA. ( 

z 

(II-2.18) 

K ( 2tt )~ 2 1 <^0 v 1 *n^> 

1 2 

(n-2.19) 

J k„i 2 ■>“, 


(II-2.20) 

■l“o. , l , ' dn o 


(II-2.21) 

° - £ n ( Re< 

n)) 

(II-2.22) 


where the a’s have been taken outside of the R integrals because they vary slowly 
and most of their contribution comes from a narrow range of internuclear separa- 
tion R= Re (0) . 

Chen used Morse anharmonic wavefunctions for the X’s and approximated the 
resonance energy (II-2.22) by 


£(m) _ £ (m) £(m) 

nA. nO T nO T 




x( n ) 


y\ 


(n-3.3) 
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where /2 is the zero point energy of the corresponding harmonic oscillator 
and the last term is an anharmonic correction. The /3’s are found from the set 
of linear equations (II-2.19). The cross sections are then determined from 
(H-2.18). 

The excited field for resonance was identified as the x 1 2” of N . Once the 

g 2 

excited state is identified the energy scale and averaged peak widths are shifted 
and altered by varying two parameters and until adjustment to 

the experimental data is achieved. The values giving the best fit are 

c ( nV + A(»> = 1.89 ev 
P>> = 0.152 e v 

These values are then used in the calculation of the cross-sections for excitation 
and de-excitation for other levels. (See reference (30) for the corrected figures 
for reference (29).) 

A1.3 Calculation of reaction rates; 

The reaction rate I. . , is obtained from the cross section, a. , by 

ij ij 

X ij = jjj V e CT ij ( V e) f ( V e) d3 V 3 (Al-18) 

where 

v e is the electron velocity 

f (v e ) is the distribution function for the electrons (assumed Maxwell- Boltzmann). 
Equation (Al-18) when expressed in terms of electron energy, E in Ev, becomes 
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E a.. (E) e -E/kTe dE (Al-19) 

where m e is the electron mass in grams. 

Te is the electron temperature in °K - 
k is Boltzmann's constant. 

Chen presents, in graphical form, the cross sections which allow calculation of 
(Al-19) for all transitions among the first five excited levels of N 2 . Because we are 
concerned with moderate electron temperatures, the low energy regions of the cross 
sections are of particular interest. By using (Al-19) and detailed balancing argu- 
ments we may relate the cross sections for excitation and de-excitation in the fol- 
lowing way 


/ 9 \ 3/2 1 

1,,= (1.26574, -M 


V <E - £ ‘j> = inri- "u (E) (A1 ' 20) 

where e . is the threshold for the j - i transition. Thus, near threshold we should 
use de-excitation cross sections in (Al-19) because the low energy tails are more 
readily observed. The principal of detailed balance allows calculation of the exci- 
tation reaction rate from 


X ii = X ij 


_(Ej-Ei)/kT e 


(Al-21) 


The cross-section graphs in Chen's paper (30) were digitized and numerical inte- 
gration of (Al-19) was performed for electron temperatures between 500°K and 
10000°K. The identification of the cross sections used are given in Table (1). The 
calculated reaction rates are plotted in Figures Al-1 through Al-15 as points. 
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REACTION RATE CONSTANT FORa 51 



ELECTRON TEMPERATURE (°K) 

Figure Al-9. Reaction rate constant calculated from 
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It was found that the reaction rates could be adequately reproduced by a four 
parameter function of the electron temperature of the form 

I n =Exp (p. (Al-22) 

\‘ T. 1* T lJ 

The four parameters for each cross section of Table (1) are given in Table (2). 
Equation (Al-22) reproduces the integrated cross section results, over the elec- 
tron temperature range, to better than 40% in all cases. This was considered ac- 
ceptable considering the order of magnitude accuracy of the basic cross section 
calculations. The reaction rates from (Al-22) are also plotted as lines in 
Figures Al-1 through Al-15. Equations (Al-22) and (Al-21) were used to calcu- 


late the reaction rates used in further calculations. 



APPENDIX II 


SIMPLIFICATION OF THE CONTINUITY EQUATIONS 

From section 2.4 we have for the continuity equation of the number of N* in 
level i 


Bn. _ v — i 

1 + — (F. ) = / {A. . n. - A. . n.) 

B t B z v *'* Z—i lJ J J 1 y 


j = o 


P l. 0:0.1 + 1} ^ + 1 “ "i n i ] 


i = o 


P i. 0:0.1 ^ (i + !) ^ j [n i+1 n JM - n. n.] 


j = l 


+ c i,o v ' n [0] [i n i-i + (*• + !) e ® /T n i+1 


+ (i e 8/r + i + 1 ) m] ■ 


( 2 - 22 ) 


These are a set of 6 coupled, non-linear, partial differential equations. The co- 
efficients are functions of altitude through the altitudinal dependence of the model 
atmosphere and ionosphere parameters. 

Each solution to equations (2-22) will be written as a sum of a Boltzmann dis- 
tributed density, n. , and a remainder term, e. , to simplify the coupling between 
the equations. Thus, 


and we require that 


n. = n. + e. 
11 1 


z>-z> 


(A2-1) 


(A2-2) 
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Z>=Z>=- 


(A2-3) 


i = 0 i = 0 


Equation (A2-2) and (A2-3), respectively, require the vibrational energy and density 
of the Maxwell- Boltzmann distributed density to be the energy and density of the 
solution for the set of equations (2-22). Thus, from equation (A2-1) we have using 
(A2-2) and (A2-3) 


z>- 


i =0 


Z>° 

i = 0 

Substituting (A2-1) into the flux term (A2-21) we obtain 


(A2-4) 


(A2-5) 


where 


F| = 

p; + 3 ! 

(A2-6) 


n.' 1 

(A2-7) 

F' = - D 

Vn.' +— k 
1 H J 


r e. .i 


9 ; = - d 

1 

J* 

+ 

■H 

vu 

>J 

(A2-8) 


Substituting (A2-1) into (2-22) and rearranging terms we have 




j = 0 


- A. (n.' + €.)} +P 1>0s0il v i (j + 1) [n.'.j n.' +1 

j *0 
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We have 


n i n j + e i-i n j' + i + e j + 1 . n i'-i " e i n l 


- e. n.' + e. , e. - e. e.l 

) 1 1-1 J + l 1 


p i. o:o.i v C 1 + ^J^ni+i n j'-i - n/ 


+ e i+i n j-i + e j-i n i+i - e i n j - e j n i 


+ e i+i e j-x " e j e jJ + C 1>0 v'n[0] [i (n!_ t + e.^) 

+ (i + 1) e* /T (n.' +1 + e. +1 ) - (i e e/T + i + 1) (n' L + e.)] (A2-9) 


n.' , n.' - n.' n.' = n ' n.' , - n'. n'. 

i-l j+l i i i+l j-l i j 


Using (A2-10) in (A2-9) and rearranging terms we have 


(A2-10) 


17 + bT <ft>. + 77 + £ &>. * E {A “ <">' + £ -> 

j “ 0 

+ A. (n' + €.)} +P 1>0!0il V i 6^2] (j + l)n; +1 

j«0 

- (i + 1) e A 22 i n / + 1 "i'-i 22 (i + e i + i " + !) n 'i 2] ’ € j 

j = i j =o j = l 

- 1 € j (j + v n j + (i + g i + 1 2"! j n j-i ■ 1 n i 22 (i + 
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+ (i + 1) n/ 


i + l 


Vi + i£ i-.Z] 


(j + 1) e 


i + i 


i = o 


- (i + 1) 6- j e. - i e. (j + 1) e. + (i + 1) e. + , j e j-i 
j = l j =o j = l 


(A2-11) 


+ c i,o v ' n ^ [i (°i-i + e i-i) + C 1 + !) e9/J ( n i+i + e i+i> 


- (i e 9/T + i + 1) (n.' + €.)] 

Using equations (A2-4) and (A2-5) in (A2-11) the 5th, 6th, 9th, 10th, 11th, 12th, 13th, 
and 14th terms on the right hand side are zero. Thus we have 


3 n.' 

I 

Tt 


1 + r, (f;> ‘ + TT + TZ ®‘>- = H <A “ (n ‘ * V 


j =0 


-A. (n.' + 6.)> +p i>0:0il V 


{i e._ 1 - (i + 1) e.} i n j 

j = 1 


+ {(i + 1) e i+ i - 1 e i> 2^ j n i-i 


j = 1 


+ C li0 v' n [0] [i (n.'_j + e.^) 


+ (i + 1) e * /T ( n i +1 + e i+1 ) " (i ^ /T + i + 1) 


(n* + e.)] (A2-12) 

A .2.1 The equation for the Boltzmann part of the solution. 


Multiply (A2-12) by i and sum over i to get 
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i = 0 'i=0 / z i = 0 j = 0 

fz>(z: 

L j = 1 i * 0 


- A ji ( n l + £ i» +P 1. 0:0.1 


[i 2 6. 


i-1 


- (i + 1) e il| + Zj j n i-l | Z] [i (i + € i + l " i2 6 i 1 jj 

+ C 10 v' n [0] f 2^ u 2 n.'.j + i (i + 1) e e/T n' +1 
L i = o 

- i (i e l9/T + i + 1) n.'} + 22 ^ i2 € i-i " i + *) € i^ 

i = o 

+ <? n £ {i + !) e i+i - iJ «.>] 

i = 0 j 


(A2-13) 


where we have used (A2-4) and 


£> 3 ,- -£*»[*«. 4^1 

i - 0 i = 0 L J 

=- D [ v Z] i6 i + E i ¥ t ] 

L i=0 i = 0 J 


= 0 


Consider 


22 [i2 e i-i - i ci + 1) ^.] = Z2 i2 e i-i - 22 i2 6i ■ H 1 €j 

i = 0 i = 1 i = 0 i = 0 

=E (i + l)! 5 ‘-L i2 e * 

i = 0 i * 0 


(A2-14) 



Ill 


= (2 i + 1) e i 


i = 0 


= 0 


Also consider 


22 [i ^ + e i+i ■ i2 e i ] = 22 1 (i " l} 6i ~ 22 i2 e * 

i = 0 


i - 0 


i = 0 


-Z> 

i = 0 


= 0 


Using (A2-15) and (A2-16) in (A2-13) we have 


i = 0 \ i = 0 



i= 0 


+ i (i + 1) e 0/r n.' +1 - i (i e d/T + i + 1) n.'} 


Consider 


{i 2 n.'_ x + i (i + 1) n' +1 — i (i e e/T + i + 1) n.'} 

i = 0 

= y ' {( i + l) 2 nf - i 2 n! - i n! + e d/r {i (i - 1) n! - i 2 n/}} 
i E 0 


(A2-15) 


(A2-16) 


(A2-17) 
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= {(i + 1) - i e*/T } n f 


(A2-18) 


Therefore equation (A2-17) becomes 


=ZD (A “ (n ; te < : 


i = 0 j = 0 


- A jA (n! + e.)} + C 1>0 v' n [0] ( i + 1) n! 

__ i = 0 


z> ; 


Let Z p be the partition function for the Boltzmann distributed part of the solution. 
Then, 


1 - e" 


and T v is the vibrational temperature. Therefore 




i B 0 P 


S (1 - e -a ) 2 


= n (Z p - 1) 


(A2-23) 
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Also 


i = 0 L i = o i = 0 J 

f n (Z - 1) 4 

= -D [Vn(Z p -l) + kj 


2: 1 22 {Aij (nj> + 6j) " Aji (n * + €i)} = (n ° + €o) |Z] 1 Ai ° (1 ■ 8io) 


i = 0 j = 0 


+ ( n i + e l) 





+ K * £ 2 > ( 2 ] 1 A ‘= (1 - 8 ‘^ - 2 Z 2 A ‘ 2 (1 ' Si2) f 

L i = 0 j =0 J 


+ + e k) J 7 iA ik (l-8 ik ) 


- 2] k *!<■ (i - *7 

j =0 J 


K + e o> 1 22 1 (l ■ §io) Ai ° + (n * + £i) 22 (i ~ 1) Aii l 

[i'O i“0 J 


+ ( n 2 + e 2 ) «< / ( i - 2) Aj 2 v + • • • 


+ K + e k) J 7 . (i - k ) A jk i + 


(n; + e.) (i - j) A. 


i = 0 j * 0 
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Using (A2-23), (A2-24) and (A2-25) in equation (A2-19) we have 


— n (Z - 1 ) 
B t p 


+ ^{- D 


3 n (Z D - 1 ) 

3 -L„CZ p -i) t _L_ 


v -1 r> 


+ e i)(i- j)A ij 


i = 0 j =0 


Let 


and use 


+ C l 0 v' n [Q] [n (Z p - 1 ) + n - e 6 ^ n (Z - 1 ) 


+ = (Z P - !) 


and 


B n 
B t 


= 0 


(Z - 1 ) |-^ = - — Zp ~ X) 
p B z H 


in equation (A2-26) to obtain 


"Tt + 3T (' d "It) = Z] Z] (n ' +£ i } (i ‘ i,A ; 


i = 0 j = 0 


+ C l>0 v' n [ 0 ] n [</r + 1 - e e/T \p] 


Using 


n.' =— (e _a y 
i Z v ' 


- n / 4, y 
4 > + i \'/' + 1/ 


(A2-26) 


(A2-27) 


(A2-28) 


(A2-29) 
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and dividing equation (A2-28) by n we have 

3 'P J_ / B D n\ B </> _ D B <// 2 _ V -1 V -1 T 1_ / V + Jb 

B t _ n y Bz J Bz 3 z 2 Z_j Z_j \4> + 1 \*Jj + lj n _ 

(i - j) A tj +C li0 v'-n [0] ' : + 1 - e 5/T </»]• (A2-30) 

The boundary conditions become for: set 1; 
initial condition: 0 (z, 0) = 0 

lower boundary: 0 (120, t) = 0 

Bn. n. Bn/ Be. n/ e 

upper boundary: — - = --7-=^-— L + -7— 1 - ~ 

Bz HBzBz HH 

Using (A2-23) 



Using (A2-23) and (A2-27) we have 


Thus, 
set 2; 

initial condition: 


U I 

r— n 
B z 


z 


8 00 



B i p 
B z 


= 0 

Z = 0 


0 (Z, 0) 


e -0/355 
1 _ e -l9 / 35 5 


(A2-31) 


= 7.329241 x 10' 5 
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lower boundary: 4> (0, t) = 7.329241 x 10 -5 

upper boundary: -Ll L = 0 (A2-32) 

B z 

It is worthy of note that the only coupling between equation (A2-30) and the other 
equations is through the e. /n term of the source of N* due to electron collisions. 
This fact, and the expectation that e. /n is small with respect to 1 for all j of interest, 
will be exploited in the solution of (A2-30). 


A2.2 The equations for the deviation parts of the solution. 

Returning to equation (A2-12) for the e. ’ s and treating the n'’s as known 
quantities we have 


3 

B t 


e. 

1 


3 


+ tl = £ {A “ (n ‘ + ‘i> - A f* K + - j r 

i = o 

^i')* +P i. 0:0.1 v pi e i-i - (i + i) e i> 22 j n i 

j s 1 


+ {(i + 1) e i+ i - 


i£ ‘ > E in '->] 

j - 1 J 


+ C 1,0 v> n [i ( n i'-i 


+ V,) + (i + 1) e d/T (n' +1 + e. + 1 ) - (i f'* + i + 1) 


(A2-33) 


Using equations (A2-23) and (A2-27) we have 


B 

B.t 


e. 

1 





3 = 0 


(Of + €.) - A., (n + £.)} - 


B n.' 
1 

Bt 


3 z + P l,0: 0,1 u ^i e i-l “ (i + i) e i } n ’A 
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+ {(i + 1) e i+1 - i «i> n (i/r + 1)] 

+ c 1>0 v' n [0] [i (n.'_ 1 + e._ x ) + (i + 1) e 9 ^ (n.' +1 + e. + 1 ) 
- (i e 9/T + i + 1) (n.' + e. )] 

Rearranging terms in (A2-34) and dividing by n.' , we have 


1 n. n. 

+ i> + Cj 0 v' n [0] (i eP /r + i + 1)] 


e . , n* 
1-1 1-1 


n i-l n i 


/ ^ P 1 , 0 : 0,1 V 1 +C 1,0 V> n [0] ' i] 


6 i +1 n i + l 


n. # n! 

i + i i 


[P 1,0:0,1 v n (i + 1 ) ( 0 + !) +C li0 V n [ 0 ]e‘ 


,d/T 


+ C 1(0 v* n [0] 


n: 


i-l 


n.' 

1 


(i + 1) e 9/T —ill - (i e 9/r + i + 1) 


n. 


Consider 


_9 
3 t 


thus 


e. 

i 

_ J_ 


€. 

» _ 1 

: n. # 

i 

n i 

at 

C n i ) a 

3 6 . 

i 

d 

6 . 

1 

l e i 

d t 

~ 

— 

4 - ■ __ 

d t 

n.' 

1 

n I n i 


n? 


(A2-34) 


(i + 1)] 

(A2-35) 


(A2-36) 
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Also 


Therefore 


Therefore 


e. i e. 

V_L= _ Ve. _ i_ Vn! 

(n{) 2 1 


n'. n! 
1 1 


Ve. .= n! V_!i +IiV n ' 


n? n? 
1 1 


L D /n'-lli lill 

n * Bz | 1 3z n j n! Bz + H J J 

.-i-ff-P nAIi. D !i{^ + l}l 

n' oz |_ oz n * n' L H JJ 


Let 


e. 


<P= =■ 


n: 


then (A2-36) and (A2-38) become, respectively, 


Be. 


3 , 


and 


+ ;p37 n ‘ 

Using (A2-39), (A2-40), (A2-41) in (A2-35) and rearranging terms we have 


B 

Bt 


B 2 TlB 1 f^n! n!')"l Bcp. 

Bz 2 J^n| Bz 1 n# [ Bz hJJ Bz 


8 ..) 

ij ' 


+ ^i, 0:0,1 vn + 1) >A + i) + C 10 v' n [0] (i e^ T + i + 1 


(A2-37) 


(A2-38) 


(A2-39) 


(A2-40) 


(A2-41) 
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■ L f 

: 22 A ^ 


(1 + <P,) e 


-(j-l) e/Ty 


- A *} 


(1 - S. .;) - 

v i j y 


n. 


Bn! 

1 

BF 



+ cPi-x e 5/T [P 1(0!0(1 , vniiA 


+ C 10 v ' n [0] i] 


+ V, 


i+l 


e" e/T v 


tt*l, 0:0,1 vn ( L + 1 )('P + X > 


+ C li0 v' n [0] e d/T (i + 1)] 


+ C. 


1.0 


n [0] [ i e 


<?/T v 


+ (i +.1) 


0/T-e/Ty 


- (ie‘ 


6/1 


+ 1 + 


»] 


Consider 


Thus, if y = t 


and if y = z 


n' By 


where 


n j =n 


(>/' + i) j+1 


J_ = _1_ Bn + /j j + l) 

n' By n By [\p \jj + 1J By 



i _? 2 i . i . r j uj) jv 

n'Bz H "j i// >// + lj Bz 


Using (A2-44) and (A2-45) in (A2-42) and collecting terms we have 


(A2-42) 


(A2-43) 


(A2-44) 


(A2-45) 
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— <p. - D — <p. - + D (--+2 /— - -i-til — 

Bt 1 ^ z 2 1 |_Bz \ H \<Jj 0 + 1 / BzJ Bz 


i + 1 ] ^0 _ ^0 fj_ _ j + 1* 1 fB D TJ / J_ 

_"\_0 0 + 1J Bt Bz (0 i/» + 1J > |bz + \H 

f j j + n Bi /<\1 D r i i + 11 B 2 0 D / B0\ 2 M 

1 0 0 + 1 J Bz /J \0 0 + 1 J 3 z 2 \ Bz/ \ ^2 


+ iiiiii + y a. . (i 

(0 + 1) 2 J Jl 


- S ij) + P 1, 0:0.1 vn { ( 2i + !) 0 + i> 


+ C 1>0 v'n[ 0 ] (i e* /T + i + l)j <p. = y |a.. (1 + <p.)e ° °* /Tv j. (1 - S. . ) 


(i_ i + 1 \ Bi // B 0 fi_ i + l \ f BD ^ f_l. 

^0 0 + ly Bt Bz \\p 0 + 1 ) \Bz + 1 H 


(— - -sr\ + 7^ a., (l _ 8..) 

\0 0+1/ BzJ Z__I J1 1J ' 10 0 + 1 J 

j*o 


B 2 0 

Bz 2 


_ d (Of Li. i 1 + 1 1 

\ BZ / 1 02 (0 + 1) 2 J 


+ P 1. 0:0.1 vn{ (2i + 1) 0 + i} 


+ Cj ^ 0 v' n [0] (i e d/r + i + 1) 


+ P i. 0:0.1 vn { ( 2i + 1) 0 + i> 


+ Cj ^ Q v' n [0] 


i e* /T » t(Ul) e^- s "v] + », m e 6/T v [p. M . v „. 0 


+ C[ , o t '* n [0] i] 


+ <P, 


i + 1 


e -0/T v 



o. 


vn(i +1) (0 + 1) 


+ C i.o n e ^ /T + !)] 


(A 2 - 46 ) 
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The boundary conditions for the e. ’s are, by assumption, identical for set 1 
and set 2 boundary conditions for 0 . For 


initial condition 

cp. (z, 0) = 0 

for all i 

lower boundary 

qp. (120, t) = 0 

for all i 

upper boundary 

Bn. n. 

1 __ 1 

"b7 " "IT 


thus, 

using (A2-45) 

Bn! Be. n. e. 

i _ i i 

Bz + Bz H H 


n! 

1 

r 1 fi i + l'l Bi/T| 

[_ H ^</» \p + lj BzJ 

Be. n. 

1 _ 1 

Bz H 


1 \p + 1J Bz Bz H 


using equation (A2-32) we obtain 

Bz , __ H 
- 2=00 

now using (A2-37) and (A2-39) we have 

, B Sn i 

n. — qp. + cp. 

1 Bz 1 1 Bz 

thus, using (A2-45) 


" <Pi H 


n! cp + qp. n! + (— + — \ = _ qp -2i- 

1 Bz 1 1 1 \ H \\p \p + l) Bz J H 
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again using (A2-32) we have * 


or 



(A2-47) 


From equation (A2-46) it is seen that all solutions are coupled through the 
electronic source term and that only those levels adjacent to the level of interest 
are coupled through vibration- vibrational or vibrational-translational energy ex- 
change terms. If 0 is assumed known, then equation (A2-46) represents a set of 6, 
coupled, linear, partial differential equations for the e. *s. 



APPENDIX m 


CONSTRUCTION OF FINITE DIFFERENCE EQUATIONS 

To numerically solve equations (A2-30) and (A2-46) we will use the finite dif- 
ference approximation of Crank and Nicholson (see reference (36)). This particular 
choice was made because, in a linear problem, the errors are of order two in both 
the time and space step, and the solutions are expected to be unconditionally stable. 

A3.1 Finite difference form, the </> equation. 

Equation (A2-30) becomes in six point finite difference form with 

M = i jj((l - 1) Az, (m - 1). At) l = 1, 2, • • •, L 

m = 1, 2, • • •, M (A3-1) 

where Az is the altitudinal space step size and At is the time step size 



and the subscript l and superscript m mean that the quantities are evaluated at the 

m 

altitude 120+ (-t-1) A z and time X (m-1) At, respectively. The term e. /n has been 
re tain ed in the source term in (A3-2) for completeness. Solving for the m+1 super- 
scripted terms in (A3-2) we have 
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h At , m+ i 


2(Az) 2 


c; + i + 


D ^ At ~ 

(Az) 2 




D n At 

^ +1 

2(Az) 2 * 1 


D ^ At + 1 ft Dn\ At 


2(Az) 2 n 4\ 3z A 2Az J 


*1 


+1 


+ ^1- 


+ At 


D^ A t' 

( A z) 2 




f D ^ At 


t 


1 



2(Az) 2 \ ^ z J'i 2Az 

j 4, 


-A 


w w w + 1 w + v 


(i-j)A 


i j 


Let 




1 + 


At 

( A z) 2 


= i + 2 H 



+ 





+ At 


Fo w L^ ++ 1 w + \ 



(i-j)A ij 


Thus (A3-3) becomes 


(A3-3) 


(A3-4) 


(A3-5) 


(A3-6) 


=d ? 


(A3-7) 
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The boundary conditions become 




set 1 

set 2 




m = 1 

0i = ° 

0j = 7.329241 E - 

05 



l = 1 

g 

II 

O 

0“ = 7.329241 E - 

05 


(A3 -8) 

It 

r 

0m _ 0m_ i _ o 

0m _ 0m_ i _ 0 


(A3-9) 


A3 .2 Finite difference form, the cp. equations. 


Let 




^a m = [|^ + D ( 

f 1 9 
- — +2 

1 |_3z ' 

t H 

■Rtn _ [7 i i + ^ 

\ 30 ; 

1 L\0 0 + 1 

/ "1 

+ D f-i + U 

i + 1 

\ H 10 

0 + 1. 


(A3-10) 


X* 


B0 30 fi i + l'IJ’BD 

j_0 0 + lj |Bz 


3z2 


L± + ILtll\ + T A.. (1.8..) 

W\0 2 OA + l ) 2 J 1J 11 

+ p, 0 . 0 ! vn{(2i + 1) 0 + i} + Cj 0 v'n [0] (ie 0/T + i + 1)1 (A3-1I) 

-U.« 


l y\ = 


J2 { A ij (1 + <Pj) e" 0 " 0 &/Tv (1 - S..) - 


L j =0 


+ p i , 0 :o.i 1/11 ^( 2i + 1) 0 + i) + C 1>0 v'n [0] |i e' 9/Tv 

+ (i + l)e^-^v} +Vi e* /T v [P 1O:Oii , ni 0 

+ C i,o v> ^ + ^i+i e- 0/Tv [P 1|0s0tl vn(i + 1) (0 + 1) 

+ C 10 v* n [0] e 0/T (i + 1)] 

- 'tm 


(A3-12) 
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Using (A3-10) through (A3-12) and expressing (A2-46) in finite difference 
form we have after some rearrangement of terms 


Let 


D ^ At ■Ui m+1 

D < 4, ( , 

^qpHl + l 
T 1 

D 

!£ At 

T • + 

2(Az) 2 

_(Az) 2 

2(Az) 

II 

1 

o 

r+ 

. + — 
1 A __ 


+ 

1 - 


't - 1 qj.m+ 1 


|_2(Az) 2 


2Az 




T At ^cpT 
1 T 1 


D -t At . t a . at_' 

2(Az') 2 1 2Az 


(Az) 2 


-e-i. 


+ ^y 1 ? At 


* A . = P l At 


U B. = 1 + 


2(Az) 2 
D^At 


(Az) 2 


= 1 + 2^A. 


(A3-13) 


(A3-14) 

(A3-15) 


= 


, 

[> 

r-h 

+ — 

w + 

ri-^^.v At] 

L2(Az) 2 

1 2Az 

1 

(AZ) ! ‘ J 


'qp“ (A3-16) 


D^At 


At 


L2(Az) 

Thus equation (A3-13) becomes 


- 

2 1 2Az 


+ ^y“> At 


_ 'ip^+l ~m+l 'Cg 'tmin+l _ ^ _1 m m+1 = 'V” 

i i i Y i i ^ i i 


(A3-17) 
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The boundary conditions are 




m = 1 

o 

II 

l-l -H 

©- 

for all i 


l = 1 

* cp™ = 0 

for all i 


l = L 

L cp^ - L_1 cp” = 0 

for all i 

(A3-18) 
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APPENDIX IV 


COMPUTER PROGRAM LISTINGS 


BOLTSOL 


C 


THIS PROGRAM CALCULATES THE 

VIBRATIONAL TEMPERATURE, AND NUMBER DENSITY 
OF VIBRATIONAL QUANTA PRODUCED BY SLOW ELECTRON 
COLLISIONS FOR VIBRATIONAL NITROGEN. 

THE SYMBOLOGY IS THAT CONTAINED IN THE ANALYTICAL 
FORMULATION OF THE PROBLEM. 

THE DIMENSIONS OF PSI WILL CHANGE WITH THE UPPER ALTITUDE 
LIMIT, ZUPPER, THE STEP SIZE, DELZ. 

TIMEX IS THE MAXIMUM TIME AND DELT IS THE TIME STEP. 

PSI SHOULD BE DIMENSIONED ACCORDING TO PSI ( I ALT , 2 ) WHERE 

MINIMUM VALUE OF I ALT IS 

I ALT = 1. + (ZUPPER - 120. ) /DELZ 
IPMOVE IS (ARRAY SIZE) KBYTES* ( NUMBER OF ARRAYS). 

ALL OTHER ARRAYS ARE DIMENSIONED FOR A 10. KM. DELZ. 

MMMTIN GOVERNS THE PRINTER LISTING AND PLOTTING FREQUENCY. 
TCHANG IS THE TIME AT WHICH THE STEP SIZE IS CHANGED 
TO DELT2 FROM DELT. 

ITFLAG IS A FLAG TO INDICATE TO THE PROGRAM THAT MORE 
THAN ONE TIME STEP IS USED. 

I F LAG IS RELATED TO ITFLAG AND IS INITIALLY SET TO ZERO. 
ICFLAG IS A FLAG TO INDICATE TO THE PROGRAM THAT A CORRECTION 
USING RESULTS FROM PROGRAM BOLTDEV IS TO BE APPLIED TO 
THE SOURCE TERM OF THIS PROGRAM. 


REAL *4 NED,NI , I NT , NZ ,NE ,MB AR , N2Z ,NZ 1 , NZ 2 , MB AR 1 , MB AR2 , NOZ 
LOG I CAL* 1 I MAGE ( 5 1 5 1 ) 

DIMENSION RAID! 8 ) ,TED( 80) ,ZD( 10) ,NED( 45) 

C ,GAMIZD(70UW(15),X(15),Y(15),ZP(15) 

C , ID1(15) ,ID2( 15) ,Z(69) ,NOZ ( 89 ) , RATEOE < 89 ) , RATEOD ( 89 ) 

C *TEI(2)»TE(89),NZ(6),NE(89),A(89) 

C »NZl(6)»NZ2(6),Pl(6,6)»P2(6,6)»P3(6,6)»P4(fc»fc) 

C ,INT(6,6),E(10)»B(89),C(89),D1(89),D2(89) 

C ,03(89), D( 89) ,DI ( 89) ,PSI ( 89,2 ) 

C *EC(89),F(89),TV(90),QUANTA(89),N2Z(89) 

C , TZZ ( 89 ) , ALPSI ( 90) ,GRID(8281 ) ,C0RR(89) 

C , ALN2Z ( 89 ) , ALNE ( 89 ) , ALQUAN ( 89 ) , DSOUR ( 89 ) , DSOURC ( 89 ) , T V 1 ( 1) 
X ,PHI(6»21),EPSIL(6,21),NB0LT(6,21), SUMC ( 6* 6 ) , G( 1 ) , N2Z1 ( 1 ) 
DATA INPUT. 


READ (5,1* END = 90000 ) RAID 

1 FORMAT (F5.1) 

WRITE (6,2 ) RA ID 

2 FORMAT ( ' 1' ,8F6.1 ) 

READ ( 5 ,3*END = 90000 ) TED 

3 FORMAT ( F6.0 ) 

WRITE (6,4) TED 

4 FORMAT ( • 0 • , 10F7.0 ) 


00000100 

00000200 

00000300 

00000400 

00000500 

000006.00 

00000700 

00000800 

00000900 

00001000 

00001100 

00001200 

00001300 

00001400 

00001500 

00001600 

00001700 

00001800 

00001900 

00002000 

00002100 

00002200 

00002300 

00002400 

00002500 

00002600 

00002700 

00002800 

00002900 

00003000 

00003100 

00003200 

00003300 

00003400 

00003500 

00003600 

00003700 

00003800 

00003900 

00004000 

00004100 

00004200 

00004300 

00004400 

00004500 
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READ. { 5 *5*END=90000) ZD 

5 FORMAT ( F 5 • 0 ) 

WRITE (6,6) ZD 

6 FORMAT ( 'O' * 10F6.0) 

READ (5*7* END = 90 000 ) NED 

7 FORMAT (Ell. 5) 

WRITE (6*8) NED 

8 FORMAT ( '0* *10E12.5) 

READ (5*10* E , ND =90000 ) GAM I ZD 

10 FORMAT (E13.7) 

WRITE (6*11) GAM I ZD 

11 FORMAT ( • 0 1 *8E14.7) 

READ ( 5*12*END = 90000) DEL Z * ZUPPER ,DELT t T I MEX * MMMT IN* ITFLAG 
X *TCHANG*DELT2, ICFLAG*RATEO 

12 FORMAT (F4.0*F5.0*F6.0*F/'.0*I4*I2?F7.0*F6.0*I2*E10.2) 

WRITE (6* 13) DEL Z * ZUPPER * DELT * TIMEX * MMMT IN, ITFLAG* TCHANG, DELT2* 
X I CFLAG* RATEO 

13 FORMAT ( 'O' *F4.0*F5.0*F6.0*F7.0* 14, I 2* F7.0, F6.0, 12, E 10, 2) 

READ ( 5 * 14 * END =90000 ) D I I *T I *N I , ENU , TOO 

14 FORMAT ( F6.4*F6.0*E13.7*F5.3*F5.0) 

WRITE (6*15) DII tTI*NI*ENU,TOO 

15 FORMAT ( • 0 • * F6 . 4 * F6 . 0 * E 1 3 . 7 * F 5 . 3 * F 5 . 0 ) 

READ (5*16* END = 90000 ) E 

16 FORMAT ( ( 1 X * F6 . 4 ) ) 

WRITE (6*17) E 

17 FORMAT ( 'O' * 1 0 ( 1X*F6.4) ) 

READ ( 5 ,18 * END =90000 ) (W(N)*X(N)*Y(N),ZP(N), 

C ID1( N) * ID2(N ) *N=1* 15 ) 

18 FORMAT ( 1X*4E14.6*3X*2U) 

WRITE (6*19) (W(N)*X(N)tY(N)*ZP(N)tIDl(.N)fID2(N)tN=l,15) 

19 FORMAT ( ‘O' , 1 X * 4E 14 . 6 * 3X , 2 I 1 ) 

READ (5*20* END=90000 ) R A I 

20 FORMAT (F5.1 ) 

WRITE (6*21) RAI 

21 FORMAT ( *0 1 *F5.1.) 

C RECONSTRUCTION OF THE INTEGRATED CHEN CROSS-SECTIONS 
C PARAMETERS. 

WRITE (6*26) 

26 FORMAT ( ' 0 1 * 2 X , 1 I D 1 * 2 X * 1 I I I 1 , 2X , ' J J J • * 8X , * P 1 ' * 9X , • P2 ' , 

C 1 3X * *P3* * 1 2 X * 1 P4 ' ) 

DO 2700 I M =1*15 
III = I D 1 ( I M ) + 1 
JJJ = I D 2 ( I M ) + 1 
P 1 ( I I I * J J J ) = W( IM) 

P2 ( III * JJJ ) = X( IM ) ^ 


00004600 

00004700 

00004800 

00004900 

00005000 

00005100 

00005200 

00005300 

00005400 

00005500 

00005600 

00005700 

00005800 

00005900 

00006000 

00006100 

00006200 

00006300 

00006400 

00006500 

00006600 

00006700 

00006800 

00006900 

00007000 

00007100 

00007200 

00007300 

00007400 

00007500 

00007600 

00007700 

00007800 

00007900 

00008000 

00008100 

00008200 

00008300 

00008400 

00008500 

00008600 

00008700 

00008800 

00008900 

00009000 



130 


P3( ni,jjj) = y ( im) 

PA ( III ,JJJ ) = ZP( Iii) 

2700 KITE (6,27) 1 01 ( IM ) , I D2 ( I h ) , I I I , J J J , p 1 (I I j , j j j ) 

C * P 2 ( III , J J J ) » P 3 ( 1 1 1 , J J J ) , P 4 ( 1 1 I , J J J ) 

27 FORMAT ( 1 0 * , 2 X , 2 I 1 , 3X , I 1 , 3X I 1 , 1 P4E 14 . 6 ) 

C TEST OF Ft A I FOK TABLE RANGE. 

DO 2100 1=1,7 
II = I 

2100 IF ( RA I • GE . R A I D( I ) . AND.KA I . LT . R A I D ( I + 1 ) ) GO TO 2300 
WRITE (6,22) 

22 FORMAT ('1','RAI NOT WITHIN TABLE RANGE') 

WRITE (6,23) RA I 

23 FORMAT ( 'O' ,F6.1 ) 

STOP 

2300 CONTINUE 


C SET DIAGONAL ELEMENTS OF INTEGRATED CROSS-SECTION MATRIX 
C TO ZERO. 

DO 2900 I I K = 1 , 6 
2900 INTI I IK, I IK ) = 0. 

I PMOVE = 89*4*1 
I FLAG = 0 
MLOW = 2 

BK = 8 . 61 703E-05 
THETA = 3380. 

HOAX = TCHANG/DELT 
MMM I N = 0 


DIFI . = ( D I I *N I ) / ( T I **0 • 75 ) 

CALL JACH2 ( Z»TOO,TZ,NZ,MBAR,RHOZ»HZ) 

CALCULATION OF TIME INVARIANT COEFFICIENTS AT EACH 
VERTICAL GRID POINT. 

LMAX = 1 + ( Z UP PER - 120.J/DEL2 
3900 WRITE ( 6,41 ) 

41 FORMAT ( '0',5X,'DIFF',10X,'A',12X,'8',12X,'C',12X, 
C 'D1'»12X,'D2',10X»'D3'» 5X,'L') 

DO 2000 L=1,LMAX 
2( L ) = 120 . + ( L -1 ) *DELZ 
TEST OF Z TO BE WITHIN TABLE RANGES. 

DO 2400 K= 1 , 9 


2400 IF ( Z( L).GE.ZD(K) . AND. Z ( L ) . LT . ZD ( K+l ) ) GO TO 2500 
WRITE (6,24) 

24 FORMAT ( '1', 'ALTITUDE IS NOT WITHIN TABLE RANGE') 
WRITE (6,25) Z ( L ) 

25 FORMAT ('O', F5.0) 

STOP 

2500 IK = K 


00009100 

00009200 

00009300 

00009400 

00009500 

00009600 

00009700 

00009800 

00009900 

00010000 

00010100 

00010200 

00010300 

00010400 

00010500 

00010600 

00010700 

00010800 

00010900 

00011000 

00011100 

00011200 

00011300 

00011400 

00011500 

00011600 

00011700 

00011800 

00011900 

00012000 

00012100 

00012200 

00012300 

00012400 

00012500 

00012600 

00012700 

00012800 

00012900 

00013000 

00013100 

00013200 

00013300 

00013400 

00013500 
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IL = K + 1 

C CALCULATION OF TE I • 

C TE I INDEX IS THE LINE NUMBER. 

C OF THE PLOTTED ELECTRON TEMPERATURE DATA. 

DO 2600 LL=1,2 

2600 TE I ( LL ) = GAM I ZD < IK + 10*<II-1) + L L - 1 ) 

C *(RAI - RAID(II)) + TED( ( I I -1 )* 10 +IK + LL-1) 

C CALCULATION OF 6ETAZI AND THE ELECTRON TEMPERATURE. 

BETAZI = ( TE I ( 2 ) - TE I ( 1 ) ) / ( ZD ( I L ) - ZD(IK)) 

TE ( L ) = BET AZ I* ( Z ( L ). - ZD C IK) ) + TE I ( 1 ) 

C CALCULATION OF ELECTRON DENSITY. 

DO 100 IT=1,44 

ZTEST1 = 2 0 • * ( IT-1 ) + 120. 

ZTEST2 = ZTEST1 + 20. 

IF ( Z( L) .GE.ZTESTl.AND.Zt L) .LE.ZTEST2) GO TO 200 
100 CONTINUE 

WRITE (6,51) 

51 FORMAT ( 'O' , 'ALTITUDE FOR NE CALCULATION IS OUTSIDE TABLE*) 
WRITE (6,52) Z ( L ) 

52 FORMAT ( *0* ,2X,F6.0) 

STOP 

200 ITT = IT 

GRADNE = ( NED( ITT + 1 ) - NED ( I T T ) ) / 20 . 

NE ( L ) = GRADNE*(Z<L) - ZTEST1) + N E D ( ITT) 

ZA = Z( L ) 

CALL JACH2( ZA , TOO , T Z , NZ , MB AR , RHOZ , HZ ) 

N2Z ( L ) = NZ(1) 

NOZ(L) = NZ ( 3 ) 

SUMNZ = NZ ( 1 ) + NZ( 3 ) 

TZZ(L) = TZ 

RATEOD(L) = RATEO^NOZ ( L ) 

RATEOE(L) = RATEOD(L)*EXP(-THETA/TZZ(L> ) 

C CALCULATION OF DIFFUSIONS AND B . COE FF I C I ENT S . 

DIFF = DIFI*( TZ**0.75 )/SUMNZ 

A ( L ) = (DIFF * DELT ) / ( 2 . * ( ( DE L Z* 1 . E + 05 ) **2 . ) ) 

B ( L ) = 1. + 2 • A ( L ) 

I F ( L.EO.l ) GO TO 2000 

CALCULATION OF DERIVATIVE OF PRODUCT OF DIFFUSION 
COEFFICIENT AND MOLECULAR NITROGEN DENSITY 
WITH RESPECT TO ALTITUDE. 

Z1 = Z ( L ) - 0.5-DELZ 
Z2 = Z(L) + 0 . 5*DEL Z 

CALL JACH2 ( Z 1 , TOO , T Z 1 , N Z 1 , MB AR 1 , RHOZ 1 , HZ 1 ) 

CALL JACH2 ( Z 2 , TOO, TZ2 , NZ2 , MB AR2 , RHO Z2 , HZ2 ) 

SUMNZ 1 = NZ 1 ( 1 ) + NZ1 ( 3 ) 


00013600 

00013700 

00013800 

00013900 

00014000 

00014100 

00014200 

00014300 

00014400 

00014500 

00014600 

00014700 

00014800 

00014900 

00015000 

00015100 

00015200 

00015300 

00015400 

00015500 

00015600 

00015700 

00015800 

00015900 

00016000 

00016100 

00016200 

00016201 

00016300 

00016400 

00016401 

00016402 

00016500 

00016600 

00016700 

00016800 

00016900 

00017000 

00017100 

00017200 

00017300 

00017400 

00017500 

00017600 

00017700 
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SUMNZ2 = NZ2 ( 1 ) + MZ2(3) 

DIFF1 = DIFI*(TZ1**0.75)/SUMNZ1 
DIFF2 = DIFI*( TZ2**0.75)/SUMNZ2 

DERIV = ( D I FF2 #NZ2 ( 1 ) - D I FF 1*NZ 1 ( 1 ) ) / ( DELZ* 1 . E + 05 ) 

C CALCULATION. OF C AND TIME INDEPENDENT PARTIAL D COEFFICIENTS. 

C ( L ) = (DERIV*DELT)/(2.*(DELZ*1.E+05)*NZ< 1 ) ) 

DHL) = ACL) + C ( L ) 

D2(L) = 1. - 2.*A(L) 

D3(L) = A ( L ) - C ( L ) 

WRITE (6,42) DIFF,A(L),B(L),C(L),D1(L),D2(L),D3(L),L 

42 FORMAT ( 1 O' , 1P7E13. 5,13) 

2000 CONTINUE 

cm = c < 2 ) 

Dl(l) = A ( I ) + C( 1 ) 

D2( 1 ) = 1. - 2.*A( 1 ) 

D3( 1 ) = A ( 1 ) - C ( 1 > 

WRITE (6,61) 

61 FORMAT ( • 1 • ,2X, 'ALT 1 ,6X, »RATE0E» ,8X, 'RATE0DS5X, 'N(O) * ) 

WRITE (6,62 ) ( Z ( I LN ) , RATEOE ( I LN ) , RAT EOD ( I LN ) , NOZ ( ILN) , ILN=1,LMAX) 

62 FORMAT ( ' 0 1 , OP F5 . 0 , 1 P3E 1 3 . 5 ) 

IF ( I FLAG. EQ. 1 ) GO TO 3950 

C CALCULATION OF EC AND F COEFFICIENTS FOR INITIAL TIME STEP. 

C LLL IS THE ALTITUDE STEP INDEX. 

EC ( 1 ) = 0 • 

F( 1) = 7 • 329241 E— 05 
WRITE (6,43) 

43 FORMAT ( 1 1 1 , 6X , • D I 1 OX , ' EC 1 , 10X , • F • , 1 OX , 'DENOM', 10X, 1 LLL 1 ) 

DO 3200 LLL=1,LMAX 

C SET DSOURC TO ZERO. 

DSOURC ( LLL ) = 0. 

C CALCULATION OF INTEGRATED CHEN CROSS-SECTION VALUES 
C AT ALTITUDE OF INTEREST. 

BKTE = BK*TE(LLL) 

RTE = 1 ./TE ( LLL ) 

RTE2 = RTE*RTE 
RTE3 = RTE2*RTE 
DO 2800 I N= 1 , 1 5 
IIJ =“ ID1( IN) + 1 
JJI = ID2( IN) + 1 

INT (IIJ, JJI) = EXP(P1( IIJ, JJI ) + P2( IIJ, JJI )*RTE 
C + P3( I I J, JJI >*RTE2 + P4( IIJ, JJI )*RTE3) 

2800 IMT(JJI,IIJ) = I NT ( I I J, JJI )*EXP< -( E( JJI ). - E ( 1 1 J ) ) /BKTE ) 

C CALCULATION OF INITIAL SOURCE TERM DI. 

SUMS I = 0. 

DO 3400 I I L = 1,6 


00017800 

00017900 

00018000 

00018100 

00018200 

00018300 

00018400 

00018500 

00018600 

00018700 

00018800 

00018900 

00019000 

00019100 

00019200 

00019300 

00019301 

00019302 

00019303 

00019304 

00019400 

00019500 

00019600 

00019700 

00019800 

00019900 

00020000 

00020100 

00020200 

00020300 
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SUM = <IIL-1)*INTUIL,1)*NE(LLL) 

3400 SUMSI = SUMSI + SUM 

DI(LLL) = SUMS I *DELT + RATEOE ( LL L ) *DE LT* ( F ( 1 ) + 1. - 
X F( 1)*EXP(THETA/TZZ(LLL>)) 

C END CALCULATION OF INITIAL SOURCE TERM AT ALTITUDE DENOTED BY LLL. 
IF (LtL.EQ.l) GO TO 3200 
DENOM = B(LLL) - A ( L LL ) *EC ( LLL- 1) 

EC ( L L L ) = A ( LL L ) /DENOM 

F ( LLL ) = ( D I ( LLL ) + A( LLL )*F( LLL -1))/DEN0M 
WRITE (6,44) DI(LLL),EC(LLL),F(LLL), DENOM, LLL 
44 FORMAT ('O' ,1P4E13*6, 13)- 
3200 CONTINUE 

C CALCULATION OF INITIAL PSI. 

MMM =' 1 

LLMAX = LMAX - 1 

PS I ( LMAX- 1,1) = F ( LMAX- 1 )/( 1 • — E C ( LMAX- 1 ) ) 

DO 3500 NNN = 2 , LLMAX 
NNNN = LMAX - NNN 

PS I ( NNNN , 1 ) = EC(NNNN)*PSI(NNNN+1,1) + F ( NNNN ) 

NNNNA = NNNN + 1 

TV ( NNNNA ) = THETA/ALOG( ( PSI ( NNNNA, 1 > + 1 . > /P S I ( NNNNA , 1 )) 

OUANT A ( NNNNA ) = N2Z ( NNNNA ) *PS I ( NNNNA , 1 ) 

3500 CONTINUE 

PSI (LMAX, I) = PS I ( LMAX - 1,1) 

PS I ( LMAX + 1 , 1 ) = PS I ( LMAX , 1 ) 

TV (LMAX ) = THETA/ALOG( ( PSI (LMAX , 1 ) + 1 . ) / PS I ( LMAX , 1 )) 
TV( 1) = 0. 

QUANTA ( 1 ) = 0. 

QUANTA ( LMAX ) = N2Z(LMAX )*PSI(LMAX , 1 ) 

C WRITE RESULTS ON OUTPUT TAPE* 

WRITE (10) ( Z( ILN) ,TZZ( ILN) ,TE( ILN) ,TV( ILN) ,N2Z( ILN) ,NE( ILN) , 

C QUANTA ( ILN) ,PSI ( ILN, 1 ) ,DI ( ILN ) , D ( I LN ) , MMM , DE L Z , 

X MMAX ,DELT ,'I FL AG, I LN = 1 , LMAX ) 

WRITE (6,34) 

34 FORMAT ( « 1 * ,2X , * ALT * , 2X , 1 T * , 6X , ■ T E 1 , 

C 4X , * T V 1 , 4X , 1 N ( N2 ) 1 , 8X,*NE*, 12X , 1 QUANTA • , 

C OX^PSI 1 , 5X , 1 SOURCE I NT * , 5X , * SOURCE * , 5X , 1 EL EC SOURCE*, 

X 2X , 1 CORRECT I ON 1 ) 

WRITE (6,35) >MMM, MMAX, DELT , IFLAG, I C F L AG , MMMI N 

35 FORMAT ( • 0 1 , 6X , 1 4 , 6X , 1 4 , 6X , F6 • 0 , 6X , I 2 , I 2 , I 4 ) 

WRITE (6,28) ( Z( IKN ) ,TZZ( IKN) ,TE( IKN) ,TV( IKN) ,N2Z( I KIM) ,NE( IKN) 

C , QUANTA ( IKN) , PSI ( IKN, 1 ) , D I ( I KN ) , D ( I KN ) ,DSOUR( IKN) , 

X DSOURC ( IKN) f IKN* 1, LMAX) 

28 FORMAT (• * , 0PF5 . 0, 0P3 F6 . 0 , 1 P8E 1 3 . 5 ) 

C CALCULATION OF COEFFICIENTS AND PSI FOR OTHER TIME STEPS. 
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C MMMT IS AN INDEX CONTROLLING THE AMOUNT OF PRINTED OUTPUT • 

C MMM IS THE TIME STEP INDEX. 

MMMT = 2 

IF ( ICFLAG.EG.O) GO TO 3950 

READ ( 1 1 * END=90001 * ERR = 90002 ) MMM IN * DE LT 1 , L BOUND, DEL Z 1 * DELZ2 * 
X LINC,(G( 1 ) f N2Zl( 1 ) tTVIC 1 ) *< I t PHI ( I , ILN) , EPSIL( 1 1 ILN) , 

X NBOLT { I * ILN) ,I = 1,6*1 ), ILN = 1* LMAX) 

3950 DO 3600 MMM=ML OW ,MM AX 
DO 3700 LLLL=1*LMAX 

IF (LLLL.EO.l) D ( 1 ) -D1 ( 1 ) *PS I ( 2 t 1 )+DI(l) 

IF (LLLL.EO.l) GO TO 3700 

C CALCULATION OF INTEGRATED CHEN CROSS-SECT I ON VALUES 
C AT ALTITUDE OF INTEREST. 

BKTE = BK*TE<LLLL) 

RTE = I./TE(LLLL) 

RTE2 = RTE*RTE 
RTE3 = RTE2*RTE 
DO 3801 I JK = 1 *15 
INI = ID1UJK) + 1 
INJ = ID2( IJK) + 1 

I NT ( I N I 9 I N J ) = EXP< PlUNIt INJ) + P2 ( I N I , I N J ) *RTE 
C + P3IINI * INJ)*RTE2 + P4 ( I N I * I N J ) *RTE3 ) 

3801 INT(INJtlNI) = I NT < INI*INJ)*EXPI-{ E( INJ) - E ( I NI ) ) /BKTE ) 

C CALCULATION OF D COEFFICIENT AT ALTITUDE LLLL. 

SUMST = 0. 

DO 3901 1111=1*6 
DO 3901 JJJJ=1*6 

SUMS = 1 •/ ( PS I ( LLLL t 1 ) + 1 . ) *( < PS I < LLLL » . 1 ‘ )/ 

C ( PSI ( LLLL* 1 )+l.))**(JJJJ-l) )*(IIII-JJJJ) 

C ^INT( III I *JJJJ )*NE(LLLL) 

3901 SUMST = SUMST + SUMS 

IF ( ICFLAG.EQ.O) GO TO 1500 
1100 I F(MMM-MMMIN) 1500*1200*90003 
1200 SUMT = 0. 

DO 1300 11=1*6 
DO 1300 J J = 1 *6 

SUMC ( 1 1 * J J ) = INTI 1 1 * JJ )*NE < LLLL )*EPSIL ( JJ»LLLL )*< I I-JJ ) 

1300 SUMT = SUMT + SUMC(II*JJ) 

CORR(LLLL) = SUMT/N2Z ( LLLL ) 

DSOURCI LLLL ) = CORR ( LLLL ) *DELT 
IF ( LLLL.EO.LMAX ) GO TO 1000 
GO TO 1600 

1000 READ ( 11* END = 90001* ERR = 90002 ) MMM IN * DELT 1 * L BOUND, DELZ 1 * DELZ2* 
X LINC*(G( 1 ) *N2Z1U) *TV1( 1 ) , ( I * PH I ( I * I LN ) , EPS I L ( I* ILN)* 

X NBOLTI I* ILN) *I=1*6*1 ) * ILN=1,LMAX) 
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GO TO 1600 

1500 DSOURC ( ILLL ) = 0. 

1600 CONTINUE 

DSOUR(LLLL) = $UMST*DELT + DSOURC ( LLLL ) 

D(LLLL) = D1(LLLL>*PSI(LLLL+1* 1 ) + D2 ( LLLL ) *PS I ( LLLL * 1 ) 

C + Q3(LLLL)*PSI(LLLL-1* 1 ) + DSOUR ( LLLL ) 

X + RATE0E(LLLL)*DELT*(PSI(LLLL,1) + 1. 

X - PSI (LLLL*1)*EXP(THETA/TZZ( LLLL ) ) ) 

DEN0M1 = B(LLLL) - A ( LLLL ) *EC.( LLLL-1 ) 

EC ( LLLL ) = A ( LLLL ) /DEN0M1 

F ( LLLL ) = ( D { LLLL ) + A<LLLL)*F( LLLL-1) ) /DEN0M1 
3700 CONTINUE 

C CALCULATION OF PSI* VIBRATIONAL TEMPER A TURE , ETC . 

PSKLMAX-1* 2 ) = F( LMAX-1 )/( l.-EC( LMAX-1 ) ) 

PSKLMAX* 2 ) = PSKLMAX-l* 2 ) 

PSI ( LMAX+1* 2 ) = PSKLMAX, 2 ) 

DO 4000 INN=2*LLMAX 
NNNNM = LMAX - INN 

PSKNNNNN* 2 ) = EC ( NNNNN ) #PS I ( NNNNN+ 1 * 2 ) + F(NNNNN) 

NNNNM = NNNNN + 1 

TV ( NNNNM ) = THETA/ALOG( ( PSKNNNNM* 2 )+ 1 .)/ PS I< NNNNM * 2 )) 

4000 QUANTA ( NNNNM ) = N2Z ( NNNNM ) *PS I( NNNNM * 2 > 

TV ( LMAX ) = THETA/ALOG( ( PSKLMAX * 2 ) +1 . ) / PS I( LMAX , 2 )) 

T V { 1 ) = 0. 

QUANTA ( 1 ) = 0. 

QUANTAiLMAX ) = N2Z(LMAX )*PSI(LMAX * 2 ) 

C WRITE RESULTS ON OUTPUT TAPE. 

WRITE (10) ( Z( ILN) *TZZ( ILN) *TE( ILN)*TV(ILN) *N2Z(ILN) *NE( ILN) , 

C QUANTA( ILN) *PSI( ILN* 2 ) * D I( I LN ) * D ( I LN ) * MMM * DEL Z * 

X MMAX, DELT * I FLAG* ILN=1*LMAX ) 

C TEST FOR PRINTING AND PLOTTING. 

I F ( MMM.NE.MMMT ) GO TO 3600 
MMMT = MMMT + MMMTIN 
WRITE (6*34) 

WRITE (6*35) M MM * MM A X , D E L T , IFLAG* ICFLAG*MMMIN 

WRITE (6*28) ( Z( IKN) *TZZ( IKN) *TE( IKN ) *TV( IKN) ,N2Z ( IKN) *NE(IKN) 

C *OUANTA( IKN) *PSI( IKN, 2 ) * D L( I KN ) * D ( I KN ) * DSOUR ( I KN ) * 

C DSOURC ( IKN) * I KN= 1 * LMAX ) 

C CONSTRUCTION OF PRINTER PLOTTER ARRAYS. 

DO 4100 I MN = 2 * L MAX 

ALPSI(IMN) = ALQG1 0 ( PS I( I MN * 2 )) 

ALN2Z ( I MN ) = AL0G10 ( N2Z ( IMN ) ) 

ALQUAN(IMN) = ALOG 1 0 ( QUANT A ( I MN ) ) 

4100 ALNE(IMN) = AL0G1 0 ( NE ( I MN ) ) 
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WRITE (6,29) 

29 FORMAT ( ' 1'., 4-6 X, 'TEMPERATURES VEkSUS ALTITUDE') 

XR = 1000. 

XL = 0. 

YT = 10000. 

YB = 0. 

NUMPR1 = LlMAX + 1 
NUMPR2 = LMAX 

CALL PL0T2 ( GR I D , XR , XL , YT , YB ) 

CALL PLOT 3 ( ' T ' , Z , T ZZ ,NUMP R 1 ) 

CALL PL0T3 ( ' V 1 , 2 , TV, NUMPR1 ) 

CALL PL0T4 ( 20 ,' TEMPERATURE DEG. K ') 

WRITE (6,30) 

30 FORMAT! 'O' ,59X, 'ALTITUDE (KM.)') 

WRITE (6,31) 

31 FORMAT ( ' 1 ' ,51X, ' PS I VERSUS ALTITUDE') 

YT = 0.0 

YB = -5.0 

CALL PL0T2 ( GR I D , XR , XL , YT , YB ) 

CALL PL0T3 ( * * ' , Z , ALPS I ,NUMPR2 ) 

CALL PL0T4 (3,'PSI' ) 

WRITE (6,30) 

YT = +9.0 
YB = +4.0 
WRITE (6,32) 

32 FORMAT ( 1 1 ', 46X QUANTA AND ELECTRON DENSITIES') 

CALL PLOT 2 ( GR I D , XR , XL , YT , YB ) 

CALL PL0T3 ( ' E ' , Z , ALNE , NUMPR2 > 

CALL PL0T3 ( ' Q ’ , Z , ALOUAN , NUMPR2 ) 

CALL PLOT 4 (23, 'NUMBER DENSITY CM**-3 ') 

WRITE (6,30) 

3600 CALL CMOVEf IPM0VE,PSI( 1,2) , 1 , P S I ( 1, 1) ,1) 

CONTINUE 

C TEST FOR CHANGE OF TIME STEP. 

MLOW = MMAX + 1 

IF ( ITFLAG.EO.l .AND. IFLAG.EO.O) GO TO 3800 

END FILE 10 

STOP 

C SET UP DO LOOP INDICIES FOR SECOND TIME INCREMENT STEP SIZE. 
3800 DELT = DELT2 

MMAX = (TIMEX - TCHANG ) /DELT2 + MLOW - 1 
I FLAG = ITFLAG 
GO TO 3900 
90000 WRITE (6,33) 

33 FORMAT ( ' 1 ' , ' END OF FILE ENCOUNTERED ON CARD READER') 
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STOP 

90001 WRITE (6,84) 

84 FORMAT ( 1 1 1 , ■ END OF FILE ENCOUNTERED ON INPUT TAPE') 
STOP 

90002 WRITE (6*85) 

85 FORMAT ('1»*'I0 ERROR ENCOUNTERED ON INPUT TAPE') 
STOP 

90003 WRITE (6*86) 

86 FORMAT {'l'*'MMMIN IS LESS THAN MMM ' ) 

STOP 

END 
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BOLTDEV 


C THIS PROGRAM CALCULATES THE DEPARTURES IN NUMBER DENSITY 
C FROM A BOLTZMAN DISTRIBUTION OF VIBRATIONAL QUANTA 

C PRODUCED BY AN ELECTRONIC SOURCE. 

C THE SYMBOLOGY IS THAT CONTAINED IN THE ANALYTICAL 
C FORMULATION OF THE PROBLEM. 

C THE DIMENSIONS OF PSI WILL CHANGE WITH THE UPPER ALTITUDE 
C LIMIT, ZUPPER, THE STEP S I ZE , DE L Z , TH E MAXIMUM TIME, TIMEX, 

C AND THE TIME STEP SIZE,DELT. 

C PHI SHOULD BE DIMENSIONED ACCORDING TO PH II 6 , I ALT ) WHERE 
C MINIMUM VALUE OF I ALT IS 

C I ALT = 1. + (ZUPPER - 120. ) / D E L Z 

C ALL OTHER ARRAYS ARE DIMENSIONED FOR A 10. KM. DELZ. 

C MMMT IS AN INDEX CONTROLLING THE AMOUNT OF PRINTED OUTPUT. 

C MMMT I N GOVERNS THE PRINTER LISTING AND PLOTTING FREQUENCY. 

C ISTART IS THE TIME STEP FROM BOLTSOL RESULTS AT WHICH THE 

C CALCULATIONS OF DEVIATIONS IS STARTED. 

C LINC CONTROLS THE ALTITUDE INCREMENT FOR THE DEVIATION 
C ' PART OF THE PROGRAM, EQUIVALENT DELZ =LINC*DELZ 
C DELZ IS THE ALTITUDE INCREMENT FOR THE BOLTZMAN PART 
C OF THE PROBLEM. 

C DELT IS THE TIME INCREMENT FOR THE BOLTZMAN PART 
C OF THE PROBLEM. 

C EDELT IS THE EQUIVALENT TIME STEP FOR THE DEVIATION 
C PART OF THE PROBLEM. 

C TINC IS THE NUMBER OF DELT INCREMENTS IN EDELT. 

C LBOUND IS THE LOWER BOUNDARY FOR THE DEVIATION CALCULATION. 

REAL *4 N E D > N I , I NT ,N Z , NE , Mb AR ,N2 Z , NZ 1 , NZ 2 , MB AR 1 , MB AR2 , NU , NBOL T 
X , NOZ 

INTEGER TINC ,T INC2 
LOG I C A L* 1 I MAGE (5151 ) 

DIMENSION Z( 89 ) ,RAID(8) ,TED< 80) , ZD ( 10) , NED (45) 

X , GAM I ZD ( 70), NZ ( 6 ) ,NE ( 89 ) , A ( 89 ) , NOZ ( 89 ) 

X ,W( 15) ,X( 15 ) ,Y( 15 ) ,ZP( 15) ,ID1( 15) 

X ,ID2(15),TEI(2),TE(89),P1(6,6) 

X ,P2( 6,6) ,P3( 6,6) ,P4( 6,6 ) , INT ( 6,6) 

X ,NZ1( 6) ,NZ2( 6) ,GRADT( 89) ,HPRESS( 89) 

X ,E< 10) ,B( 89) , MM AX 1 ( 3 ) , DE LT 1 ( 3 ) , I FLAG ( 3 ) , 

X N 2Z ( 89 ) , MMM(3) ,DI (89) , P S I ( 30,3) ,RATE0E(89) ,RATE0D(89) 

X , T Z Z ( 8 9 ) ,GRID( 8281 ) ,SUME1A( 89) , SUME2 A ( 89 ) 

X ,DERIVD( 89) , D I F F ( 89 ) , H { 89 ) , DP S I Z ( 89 ) ,DPSIT(89) 

X ,D2PSIZ( 89) , ALPHA ( 6,89) , BET A (6,89 ) ,GAMMA(6,89) 

X , NU ( 8 9 ) , P 1001 (89), TV (30, 3) ,QUANTA(30,1 ) ,D(30, 1 ) 

X ,PHI(6,89),C1(6,89),C2(6,89),C3(6,89),C(6,89),F(6,89) 

X ,EC( 6,89) ,ALPHI1(89),ALPHI2( 89) ,ALPHI3(89) ,ALPHI4( 89) , 

X A L PH I 5 ( 89) ,ALPHI6( 89 ) , EPS I L ( 6 , 89 ) , NBOLT ( 6 , 89 ) , REPS I L ( 6 , 89 ) 
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C DATA INPUT. 

READ ( 5,1 ,END=90000) RAID. 

1 FORMAT (F5.1) 

WRITE ( 6 , 2 ) R A I D 

2 FORMAT ( ' 1 ' * 8 F 6 * 1 ) 

READ ( 5,3,END=90000) TED 

3 FORMAT ( F 6 • 0 ) 

WRITE (6,4) TED 

4 FORMAT ( , O l ,10F7.0) 

READ ( 5 ,5 ,END = 90000 ) ZD 

5 FORMAT (F5.0) 

WRITE ( 6 , 6 ) ZD 

6 FORMAT ( 'O' ,10F5.0> 

READ ( 5 ,7,END=90000 ) NED 

7 FORMAT (Ell. 5) 

WRITE ( 6 , 8 ) NED 

8 FORMAT ( 'O' ,10E12.5) 

READ ( 5,10,END = 90000) GAM I ZD 

10 FORMAT (E13.7) 

WRITE (6,11) GAM I ZD 

11 FORMAT ( 'O' , 8E 1 4 . 7 ) 

READ ( 5,12,END=90000) DEL Z , ZUPP ER , D E LT , T I MEX , MMMT I N , L I NC , T I NC 
X ,LB0UND,TINC2,RATE0 

12 FORMAT ( F4.0,F5.0,F6.0,F7.0, 14 ,41 3 t E 10.2 ) 

WRITE ( 6, 13 )DELZ,ZUPPER,DELT, TIMEX, MMMT IN, LING, TINC 
X , L BOUND ,T INC 2 ,RATEO 

13 FORMAT ( *0' , F4 . 0 , F5 . 0 , F 6 . 0 , F 7 . 0 , I 4 , 4 I 3 , E 1 0 . 2 ) 

READ ( 5,14,END=90000) D I I , T I , N I , ENU , TOO 

14 FORMAT (F6.4,F6.0,E13.7,F5.3,F5.0) 

WRITE (6,15) DI I ,TI ,NI , ENU, TOO 

15 FORMAT ( '0' , F6 . 4 , F6 . 0 , E 1 3 . 7 , F 5 . 3 , F 5 . 0 ) 

READ ( 5,16,END=90000) E 

16 FORMAT ( ( 1 X 7 F 6 . 4 ) ) 

WRITE (6,17) E 

17 FORMAT ( • 0 1 ,10( 1X,F6.4) ) 

READ ( 5 ,18 ,END=90000 ) ( W ( N ) , X ( N ) , Y ( N ) , ZP ( N ) , 

C ID1(N),ID2(N),N=1,15) 

18 FORMAT ( 1X,4E14.6,3X,2I1) 

WRITE (6,19) (W(N),X(N),Y(N),ZP(N),ID1(N),ID2(N),N=1,15) 

19 FORMAT ( *0* , 1 X , 4E 14 . 6 , 3 X , 2 I 1 ) 

READ ( 5,20,END = 90000) R A I 

20 FORMAT ( F5.1 ) 

WRITE ( 6,21 ) „RAI 

21 FORMAT ( 'O 1 ,F5.1 ) 

READ (5,22) ISTART 
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22 FORMAT (14) 

WRITE (6,23) ISTART 

23 FORMAT ( 'O' ,14) 

CALL JACH2I Z , TOO , T Z , NZ , MBAR , RHQZ , HZ ) 

C RECONSTRUCTION OF THE INTEGRATED CHEN CROSS-SECTIONS 
C PARAMETERS. 

WRITE (6,26) 

26 FORMAT ( ' 0 ' , 2X , • I D ' , 2X , ' I 1 1 ' , 2X , ' J J J ' , 8X , ' P 1 • , 9X , • P2 ' , 

C 13X, 'P3 ' ,12X, 1 P4 • ) 

DO 2700 I M= 1 , 1 5 
III = ID1( IM) +.1 
JJJ = ID2( IM) + 1 
P 1 ( 1 1 I , JJJ ) = W( IM) 

P2( I I I , JJJ ) = X( IM) 

P3 ( III ,JJJ ) = Y ( I M » 

P4( I I I , JJJ ) = ZP( IM ) 

2700 WRITE (6,27) I D 1 ( I K ) , I D2 ( I M ) , I I I , J J J , P 1 ( I I I , J J j ) , 

X P 2 ( III,JJJ),P3( III, JJJ) ,P4( III, JJJ) 

27 FORMAT ( ' 0 ' , 2X , 2 1 1 , 2 ( 3X , I 1 ) , 1 P4E 14. 6 ) 

DIF I =<DII*NI)/(TI**0.75) 

C IPMOV AND I MMOV ARE (ARRAY SIZE (KBYTES#! NUMBER OF ARRAYS) 
IPMOV = 30*4*2 
I MMOV = 1*4*2 
BK = 8 . 61 703E-05 
THETA = 3380. 

ISTAR = ISTART - 1 
LMA X = 1. + ( ZUPPER - 120. ) /DELZ 
C SET PHI ( I , L ) , EC ( I ,LBOUND) , AND F(I,LBOUND) TO ZERO. 

C SET REPS I L ( I , LBOUND ) AND DIAGONAL ELEMENTS OF THE 
C INTEGRATED CHEN CROSS-SECTION MATRIX TO ZERO. 

DO 3000 IJ = 1,6 
DO 3001 K=l, LBOUND 
EC ( I J ,K ) = 0. 

3001 F ( I J , K ) = 0. 

REPSILt IJ, LBOUND) = 0. 

INTI I J, I J ) = 0. 

DU 3000 LL = 1 , LMAX » L I NC 
3000 PHI ( I J , L L ) = 0. 

C READ TAPE CONTAINING BOLTZMAN PART SOLUTION. 

READ ( 10,END=90001»ERR=90002) ( Z ( I LN ) , TZZ ( I LN ) , TE ( I LN ) , 

X TV( ILN,1 ) ,N2Z( ILN) ,NE( ILN) , QUANTA ( ILN, 1 ) , 

X PS I ( ILN, 1 ) ,DI ( ILN) ,D( ILN, 1 ) ,MMM( 1 ) , DELZ1 , 

X MMAXK 1),DELT1(1),IFLA6( 1) , ILN = 1,LMAX ) 

IF (ISTAR - MMM ( 1 ) ) 90003,3600,4100 
3600 DO 3700 H = 2,3 
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00006509 
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00007900 

00007901 
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00008600 

00008700 

00008800 
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3700 READ ( 10, END=90001 , ERR=90002 ) ( Z ( I LN ) , TZZ C I LN ) , TE ( I LN ) , 

X TVULN,M),N2Z( ILN) ,NE( I LN ) , QUANTA ( I LN, 1 ) t 

X PSI ( ILN,M) ,DI < ILN) ,D< I LN , 1 ) , MMM ( M ) , DE LZ 1 , 

X MMAXKM.) ,DELT1(M) ,IFLAG(M) ,ILN=1,LMAX) 

MMMT = ISTAR + 1 
GO TO 6000 

4100 DO 1000 JK = 3 , I STAR 

1000 READ ( 10,END=90001,ERR=90002) 

MMMT = ISTART + 1 
DO 4200 LM=1,3 

4200 READ ( 10 , END = 90001 , ERR = 90002 ) ( Z ( I LN ) , TZZ < I LN ) , TE U LN ) , 

X TV(ILN,LM) ,N2Z{ILN.) ,NE< ILN ) * QUANTA ( ILN, 1 ) , 

X PSI( ILN,LM),DI< ILN),D< I LN , 1 ) , MMM ( LM ) , DE LZ 1 , 

X MMAXl(LM) ,DELT1(LM) , I FLAG ( LM ) , I LN=1 , LMAX ) 

GO TO 6000 

3900 DO 3901 0=1 ,T I NC 

CALL CMOVE ( I PMOV , PS I ( 1 , 2 ) , 1 , PS I ( 1 , 1 ) , 1 ) 

CALL CMOVE (I PMOV ,TV < 1 ,2 ) , 1 ,TV ( 1 , 1 ) , 1 ) 

CALL CMOVE ( I MMOV , MMM ( 2 ) , 1 * MMM <1 ) , 1 ) 

CALL CMOVE ( I MMOV ,MM AX 1 ( 2 ) , 1 ,MM AX 1 ( 1 ) , 1 ) 

CALL CMOVE ( I MMOV , DELT 1 ( 2 ) , 1 , DELT 1 ( 1 ) , 1 ) 

CALL CMOVE( I MMOV , I FL AG ( 2 ) , 1 , I F LAG (1) ,1 ) 

READ ( 10,END=90001,ERR=90002) ( Z ( I LN > , TZZ ( I LN ) , TE ( I LN ) , 

X TV ( ILN, 3 > ,N2Z< ILN) ,NE ( I LN ) , QUANTA ( I LN, 1 ) , 

X PSI( ILN,3),DI(ILN),D( I LN , 1 ) , MMM < 3 ) , DELZ 1 , 

X MMAX1 (3) ,DELT1(3 ), I FLAG ( 3' ) , ILN = 1,LMAX> 

IF (I FLAG { 3) •EQ.l.AND. IFLAG ( 2 ) • EQ* 0 ) GO TO 6000 

3901 CONTINUE 
GO TO 3800 

C CALCULATION OF A, 8 AND OTHER TIME INVARIANT TERMS. 

6000 DELT = DELTK3) 

IF ( IFLAG(3).E0.1), TINC = TINC2 
EDELT = DELT*TINC 
DO 2500 11=1, LMAX, LINC 

CALL JACH2 AH I I ) , TOO, T Z ,NZ , MBAR , RHOZ , HZ ) 

HPRESS ( II) = HZ ^M BAR/2 8. 

NOZ (II) = NZ { 3 ) 

SUMNZ *= NZ ( 1 ) + NZ ( 3 ) 

DIFF(II) = DI F I # ( TZ**0 .75 )/SUMNZ 

A( II). = {DIFF(II )*EDELT )/( 2.*( ( DELZ*LINC*l.E+05)^*2. ) ) 

B( I I ) = 1. + 2*A( II ) 

T05 = T Z^^O . 5 

NU (II) = 1.48209E-11*T05 

PlOOl(II) = 2.6003E-06*TZ*EXP(91.5/TZ) 

RATEOD (II) = RATEO^NOZ (II) 
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RATEOE (II) = RATEOD( II )*EXP(-THETA/TZZ("fI ) ) 

C CALCULATION OF DERIVATIVES OF D I FF . , TEMP . AND DENSITY SCALE HEIGHT, 
IF (I I .E0.1 ) GO, TO 2500 
Z1 = Z ( I I ) -0 • 5*D EL Z 
Z2 = Z ( 1 1 ) + 0.5*DELZ 

CALL JACH2 ( Z1 ,T00,TZ1 ,NZ1 ,MBAR1 ,RHOZ 1 , HZ 1 ) 

CALL JACH2 ( Z2 tTOO ,TZ2 tNZ2 ,MBAR2 »RH0Z2 ,HZ2 ) 

GRADT ( II ) = (TZ2 - TZD/DELZ 

Hill) = HPRESS ( I I ) / ( 1 + HPRESS (II) * GRADT < I I ) /TZZ (I I ) ) 

SUMNZ1 = NZ 1 ( 1 ) + NZ1 ( 3 ) 

SUMNZ2 = NZ2( 1 ) + NZ2(3> 

DIFFI = DIFI*(TZ1**0.75)/SUMNZ1 
DIFF2 = D I F I * ( TZ2 sSes!t 0.75) /SUMNZ2 
DERIVD (II) = (DIFF2 - DIFFI )/(DELZ*l.E+05) 

2500 CONTINUE 

GRADT ( 1) = GRADT ( 2 ) 

H ( I ) = HPRESS ( 1 ) / ( 1 + HPRESS(l) * GRADT ( 1 ) /TZZ ( 1 ) ) 

DER I VD ( 1 ) = DERIVDI 2 ) 

LINCP1 = LINC + LBOUND 
LLMAX = LMAX - LINC 
MMMAX = MMAX1 ( 3 ) - I 
3800 DO 4000 L = L I NCP 1 , L LMAX , L I NC 
C CALCULATION OF DERIVATIVES OF PSI WITH RESPECT TO 
C ALTITUDE AND TIME. 

DPSIZ(L) = (PSIIL + 1,2) - PSI<L-1,2) >/<2.*<DELZ*l.E+05)) 

D2PS I Z { L ) = ( PS I ( L + l t 2 ) -2.*PSI(L,2) + PSI<L-1,2))/ 

X ( { DELZ*l.E+05)**2 ) 

DELTS = DELTK1) + DELTK3) 

DPSIT(L) = ( PS I ( L 1 3 ) - PSI(L» 1 ) ) / ( 2. *DELT S > 

C CALCULATION OF INTEGRATED CHEN CROSS-SECTION VALUES 
C AT ALTITUDE OF INTEREST. 

BKTE = BK*TE(L) 

RTE = l./TE(L) 

RTE2 = RTE*RTE 
RTE3 = RTE2*RTE 
DO 2800 I N= 1 , 1 5 
IIJ = I D 1 ( IN) + 1 
JJI = I D2 ( IN) + 1 

INT (IIJ, JJI) = EXP( Pl( 1 1 J, JJI ) + P2( 1 1 J, JJI )*RTE 
C + P3( I I J,JJI )*RTE2 + P4( IIJ, JJI )*RTE3) 

2800 I NT (JJI, IIJ) = INT( I I J,JJI )*EXP( -( E( JJI ) - E ( I I J ) ) /BKTE ) 

DO 4000 1=1,6 

BRAKET = (I-1)/PSI(L,2) - I/{PSI(L,2) + 1.) 

ALPHA ( I , L ) = DERIVD(L) + D I FF < L ) * ( -1 . / < H ( L ) *1 . E+05 ) + 2. *BRAKET 
X *DPSIZ(L)) 
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C CALCULATION OF SOURCE TERMS. 

C CALCULATION OF SUMMATIONS FOR USE IN BETA AND GAMMA. 

SUMA1 =0. 

SUMA2 =0. • 

DO 5000 I I L = 1 1 6 

SUM 1 = I NT ( I tIIL)*NE(L> 

SUM A 1 = SUM A 1 + I NT ( I I L ♦ I ) *NE { L ) 

SUM 2 = SUM1*(1. + PHIt IILtL) )*EXP( (I - I I L ) *THETA/TV ( L , 2 ) ) 

5000 SUMA2 = SUMA2 + SUM 2 

BETA ( I f L ) = BRAKET*DPSIT(L) - DP S I Z ( L ) *BRAKET* ( DER I VD ( L ) + 

X DIFF(L)*(-l - ./(H(L)*l.E + 05) + BRAKET*DPS I Z ( L ) ) ) 

X - DIFF(L)*BRAKET*D2PSIZ(L> - DIFFt L >*( DPSIZtL)**2 )*<-<I-l)/ 

X (PSI(Lt2)**2 ) + I/( ( PS I ( L 1 2 ) + -l.)**2 )) 

X + SUMA1 + P1001(L)*NU(L)*N2Z(L)*(<2.*I - l.)*PSI(Lt2> + I -1.: 
X + RATEOE(L)*( ( I-1)*EXP(THETA/TZZ( L) ) + ( 1-1 ) + X.) 

IF (I .EQ.l ) GO TO 7000 
IF ( I . E 0 . 6 ) GO TO 8000 

GAMMA ( I t L ) = SUMA2 - BET A ( I * L ) + P 1001 ( L )*NU ( L >*N2Z ( L ) *< ( 

X 2.*I - l.)*PSI<L»2> + I - 1 • + < 1-1. )*PSI(L»2)*PHH I-1,L) 

X *EXP{ THETA/TV(L*2 ) > + I*(PSI(Lt2) + 1 . ) *PH I ( I .+ 1 , L ) 

X *EXP(-THETA/TV(L*2) ) ) 

X + RATEOEIL )*( ( 1-1 )*EXP(THETA/TV(L»2) >*(PHI ( I-ltL) + 1.) + 

X I*EXP(THETA/TZZ(L) - THET A/T V ( L , 2 ) ) * ( PH I ( I + 1 f L ) + 1.)) 

GO TO 4500 

7000 GAMMA(ItL) = SUMA2 - BE T A ( I * L ) + P1001 ( L )*NU ( L ) *N2Z < L ) * ( ( 

X 2.*I - 1. )*PSI(Lf2) + I - 1. + I*(PSI(L»2> + l.)*PHI(2,L> 

X *EXP(-THETA/TV(L,2) ) ) 

X + RATEOE( L )*EXP( THETA/TZZt L) - THE T A/T V ( L , 2 ) ) * ( PH I ( 2 , L ) + 1.). 
GO TO 4500 

8000 GAMMA ( I * L ) = SUMA2 - BET A ( I T L ) + P 100 1 ( L ) *NU ( L ) *N2Z { L ) * ( ( 

2 X 2.*I - l.)*PSI(Lf2) + I - 1. + (I - 1 . )*PSI(L»2)*PHI(I-1,L) 

X *EXP(THETA/TV(Lt2) ) ) 

X + RATEOE(L>*( I-1)*EXP(THETA/TV< Lf 2) )*<PHI ( I-lt L) + 1.) 

4500 CONTINUE 

C1(I,L) = A ( L ) +ALPHA( IfL)*EDELT/(2.*DELZ*LINC*l.E + 05) 

C 2 ( I t L ) = 1. - 2.*A(L> -BET A ( I * L ) *E DE LT 
C 3 ( I t L ) = A ( L ) -ALPHA! I> L ) *EDELT/.< 2. *DELZ*L INC*1# E+05 ) 

C(IfL) = Cl ( I tL ) -PH I ( I ,L + LINC)+C2( I f L )*PHI (I ,L)+ C3(I,L> 

X *PHI( I 9 L-L I NC ) +GAMM A ( I,L)*EOELT 
C CALCULATION OF EC(I»L) AMD F(I,L) COEFFICIENTS. 

DENOM = B ( L ) - A(L)*EC(IfL-LINC> 

EC ( I t L) = At L ) /DENOM 

F(ItL) = ( C ( I t L ) + A(L )*F( I,L-LINC) )/DENOM 
4000 CONTINUE 
C CALCULATION OF PHI. 
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IF (LBOUND.EG.l) GO TO 4300 
DO 4301 I L = 1 * 6 

NBOLTt ILfLBQUND) = N2Z( L BOUND ) *( 1 •/ ( PSI ( L BOUND *2 ) +!.)>* 
X ( ( PSI ( L BOUND * 2 ) / ( P S 1 1 L B OUND * 2 ) + 1 . ) ) *# tl L - 1)) 

4301 EPSIK IL»LBOUND) = 0. 


GO TO 4400 



4300 NBOLTt 1*1 ) 

= 

N 2 Z ( 1 ) 

NBOLT I 2 * 1 ) 

= 

0. 

NBOLTt 3 ,1 ) 

= 

0. 

NBOLT ( 4* 1 ) 

= 

0. 

NBOLTt 5 ,1 ) 

= 

0. 

NBOLTt 6* 1 ) 

= 

0. 

EPSILt 1 *1 ) 

= 

0. 

EPSILt 2*1 ) 

= 

0. 

EPSILt3*l ) 

= 

0. 

EPSILt4*l) 

= 

0. 

EPSILt 5*1 ) 

= 

0. 

EPSIL<6,1 ) 

= 

0. 


4400 DO 4600 11=1*6 

PHI ( I I * UMAX ) = F ( I I * UMAX )/(l. - EC(II f LLMAX )) 

PHI ( I I tlKiAX) = PHI { I I * LLMAX ) 

NBOLT( IlfLLMAX) = N2Z ( LLMAX ) * ( 1 . / ( PS I ( LLMAX , 2 ) + 1.)) 

X *( (PSI(LLMAX,2)/<PSI(LLMAX,2> + M)^(II-l)) 

NBOLT ( II* LMAX ) = N2Z{ LMAX ) * ( 1 • / ( PS I ( LMAX*2> + 1 • ) ) 

X *((PSI(LMAXt2)/(PSnLMAX f 2) + l.))**(II - D) 

EPS I L { IltLLHAX ) = PHI(II*LLMAX ) *NBOLT ( I I * LLMAX ) 

EPSIK II tLMAX ) = PHI ( II v LMAX )*NBOLT ( I I » LMAX ) 

LL LMAX = LMAX - 2*LINC 
DO 4700 L L = L I NC P 1 * LL LMAX * L I NC 
NN = LMAX - L I NC + L BOUND - LL 

PHHIIfNN) = EC( 1 1 *NN )*PHI( I I ,NN + LINC > + F( 1 1 ,NN) 

NBOLT ( 1 1 * NN ) = N2Z<NN ) * ( 1 . / ( PS I ( NN *2) + 1.)) 

X *((PSI(NN * 2 ) / ( PS I ( NN ,2) + 1 . ) ) ** < I I -1 ) ) 

4700 EPSILUIfNN > = PHI(II*NN‘ ) UNBOLT (I I * NN ) 

4600 CONTINUE 
C WRITE RESULTS TO TAPE 

WRITE ( 11 ) MMM(3) tDELTK 3) tLBOUNDf DELZt DELZl,LlNCt ( Z( ILN) t 
X N2Z( ILN),TV< ILN*2) ,< IiPHK I, ILN) » EPSIK I, ILN) f 
X NBOLT ( I * ILN) » 1 = 1 * 6* 1 ) 1 1 LN = L BOUND , LMAX * L I NC ) 

IF (MMM(2).NE.MMMT) GO TO 3500 
DO 5200 NN = L BOUND tLMAX * L INC 
SUME1 = 0, 

SOME 2 = 0. 
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nn 5100 i k = 1 1 6 

SOME = EPS I L ( IKtNN) 
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SUMEI = EPSIL(IK,NN)*(IK-1) 

SUME1 = SUMEI + SUME 
5100 SUME2 = SUME 2 + SUMEI 
SUMEI A ( NN ) = SUMEI 
SUME2A ( NN ) = SUME2 
IF ( NN • EQ * LBOUND ) GO TO 5200 
DO 5200 I K= 1 , 6 

REPSILt IK,NN ) = SUME1A(NN)/EPSIL( IK,NN) 

5200 CONTINUE 

MMMT = MMMT + MMMTIN 

91 FORMAT ( 1 1 * ,2X, *C1 < I ,L ) S2X, *C2< 1 , L ) 1 ,2X, *C3( I , L ) 1 , 5X, • C ( I ,L ) 1 , 

X 2X , 'GAMMA! I,L)S IX? *ALPHA(I»L) * , 1 X , * B ET A ( ItL) 1 , 

X 4X , 1 A ( L ) * , 4X , * B ( L ) 1 t 5X» 1 EC ( I * L ) ' » 4X , 1 F ( I, L ) * ) 

DO 4900 I JK = 1 , 6 
WRITE (6,91) 

WRITE ( 6,92) ( C 1 ( I JK,L ) ,C2( I JK,L) ,C3( I JK, L) ,C( I JK , L ) , GAMMA < I JK ,L) 
X ALPHA { I JK , L ) , BETA ( IJK,L),A(L),B(L),EC(IJK,L),F( I JK, L ) , 

X L=LINCP1,LLMAX,LINC) 

92 FORMAT { ' S1P11E10.3) 

WRITE ( 6,93) 

93 FORMAT { 1 0 1 ♦ 2X , 1 DPS I T 1 , 5X ♦ ' DER I VD 1 , 4X , 1 DPS I Z ' » 5X , 'D2PSIZ • , 8X 
X t 1 DIFF • ,5Xt 'H 1 , 8X f 'P1001 • »4X ? 'NU' ) 

WRITE (6,94) (DPSIT (L ),DERIVD(L ) ,DPSIZ(L) ,D2PSIZ(L),DIFF( L) ,H(L) 
X tP 100 1 { L ) ,NU(L)fL=LINCPltLLMAX,LINC) 

94 FORMAT ( 1 ' tlP8E10.3) 

WRITE (6,37) 

37 FORMAT ( 1 0 * ,2X, 1 MMM 1 ,1X, 1 LEVEL 1 ,3X, 1 OE L T 1 ,3X, 'EDELT* ) 

WRITE (6,35) MMM(3) , IJK, DELT, EDELT 

35 FORMAT ( ' 0 • , 2X , I 4 , 3X , I 3 , 3X , F5 . 0 , 3X , F 5 . 0 ) 

WRITE (6,34) 

34 FORMAT ( 1 0 1 , 2X , » ALT • , 2X , 1 T 1 , 4X,'TV', 

X 9X , ' N ( N2 ) 1 , 8X, » P S I 1 ,11X, 'PHI ' ,8X, 'SUM EPSILON* , 

X 3 X , 1 SUM I*EPSILON 1 ,3X, * EPSILON* ,4X, 1 BOLT DEN ' , 3X , • RE P S I L ' ) 

4900 WRITE (6,36) ( Z ( IKN ) , T ZZ ( IKN ) , TV(IKN,2), 

X N2Z( IKN) , PS I (IKN, 2 ) ,PHI ( I JK, IKN ) , SUME 1 A ( I KN ) , 

X SUME2 A ( IKN) f EPSIL( I JK f IKN ) ,NBOLT( I JK, IKN) ,KEPSIL( I JK, IKN) , 

X I KN=LBOUND, LMAX , L I NC ) 

36 FORMAT (* » ,0P F5 . 0 , 0P2 F6 . 0 , IP 5 E 14 .6 , 1P3E 1 2 . 3 ) 

C CONSTRUCITGN OF PRINTER PLOTTER ARRAYS « 

DO 4800 JJ -LBOUND , LMAX 
ALPHIl(JJ) = PHI(lfJJ) 

A L PHI 2 ( JJ ) = PHI (2, JJ ) 

ALPHI3( JJ ) = PHI ( 3, JJ ) 

ALPHI4( JJ) = PHI ( 4 , J J ) 

AL PH I 5 ( J J ) = PHI ( 5 , J J ) 
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4800 

ALPHI61JJ) = PHI ( 6 * J J ) 




00029000 

C PRINTER PLOTTING. 




00029100 


WRITE (6,29) 




00029200 

29 

FORMAT ( 1 1' * 43 Xt' DEVIATION VERSUS ALTITUDE') 


00029300 


XR = 1000. 




00029400 


XL = 0. 




00029500 


YT = +14. 




00029600 


YB = -6. 




00029700 


NUMPR1 = LMAX + 1 




00029800 


CALL PL0T2(GRID,XR,XL,YT,YB) 




00029900 


CALL PL0T31 •1 I ,Z,ALPHI1,NUMPR1) 




00030000 


CALL PL0T3( 1 2 1 , Z , AL PH I 2 , NUMPR 1 ) 




00030100 


CALL PL0T3( 1 3 1 , Z , ALPHI 3 ,NUMPR 1 ) 




00030200 


CALL PL0T3( 1 4 1 , Z , AL PH 1 4 , NUMPR 1 ) 




00030300 


CALL PL0T3 ( ' 5 • , Z , AL PH I 5 , NUMPR 1 ) 




00030400 


CALL PL0T3 ( 1 6 • , Z , ALPH I 6, NUMPR1 ) 




00030500 


CALL PLOT 4 ( 3 , 1 PH I ■ ) 




00030600 


WRITE (6,30) 




00030700 

30 

FORMAT ( *0* ,59X, 'ALTITUDE ( KM . ) * ) 




00030800 

3500 

IF ( MM M ( 2 ) ♦ L T • MMM A X ) GO TO 3900 




00030900 


END FILE 11 




00031000 


STOP 




00031100 

90000 

WRITE (6,33) 




00031200 

33 

FORMAT ( 1 1 1 , 1 END OF FILE ENCOUNTERED 

ON 

CARD 

READER 1 ) 

00031300 


STOP 




00031400 

90001 

WRITE (6,84) 




00031500 

84 

FORMAT ( 1 1 1 , 1 END OF FILE ENCOUNTERED 

ON 

INPUT 

TAPE* ) 

00031600 


STOP 




00031700 

90002 

WRITE (6,85) 




00031800 

85 

FORMAT ('l', 1 10 ERROR ENCOUNTERED ON 

INPUT TAPE ' ) 

00031900 


STOP 




00032000 

90003 

WRITE (6,86) 




00032100 

86 

FORMAT (*1','ISTART IS LESS THAN OR 

EQUAL TO 

THE FIRST TIME 

00032200 


X STEP AVAILABLE 1 ) 




00032300 


STOP 




00032400 


END 




00032500 


0396 ( 
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REACTRAT 


C 


this program calculates the reaction rates for ion-atom 

INTERCHANGE BETWEEN 0 PLUS AND N2* 

THE POPULATION OF THE V I BRAT I ON ALL Y EXCITED LEVELS 
OF N2 AND THE VIBRATIONAL TEMPERATURE COME FROM 
THE OUTPUT TAPE OF PROGRAM BOLTDEV. 

REAL * 4 Nit N2Z, NBOLT , KBOLT , KE FF ♦ KV I BL , NLTOT , KPR I M 

LOGICAL * I IMAGE! 5151 ) 

DIMENSION Z < 89 ) , N2Z(89),PHI(6,89),TV(89),EPSIL<6,89) 

X f NBOLT ( 6,89 ) , KBOLT (89), KE FF ( 89 ) , R AT I OK ( 2 I ) , K V I BL ( 6 ) 

X T I LEVEL! 6) ,TAUN ( 89 ) , ALKBOL { 2 1 ) , ALKE FF ( 2 1 ) , I I ( 6 ) , NLTOT ( 6 , 89 ) 
X , T AUB ( 8 9 ) ,TAUE( 89 ) , ALT AUB ( 2 1 > t ALTAUE < 2 1 > ♦ ALT AUN ( 2 1 ) 

X , ARRAY! 21 ) ,ARRAY1 (21 ), GRID (8281 ) 

DATA INPUT 

READ (8? 12, END= 90000 ) DE L Z , ZUPPER , DE LT , T I MEX , MMMT I N , L I NC ,.T I NC 
X , LBOUND, TINC2 

12 FORMAT ( F4.0,F5.0,F6.0,F7«0, 14,413 ) 

WRITE ( 6,13) DEL Z, ZUPPER, DELT, TIMEX, MMMT IN, L INC, TINC 
X , LBOUND, T I NC2 

13 FORMAT ( *0' ,F4.0, F5 .0 , F6 .0 , F7 . 0,14, 413) 

READ ( 5, 14, END = 90000 ) D I I , T I , N I , ENU , TOO 

14 FORMAT ( F6 • 4 , F 6 • 0 , E 1 3 • 7 , F 5 • 3 , F 5 • 0 ) 

WRITE (6,15) DI I ,TI ,NI , ENU, TOO 

15 FORMAT ( • 0 • , F6 .4 , F6 . 0 , E 13 . 7 , F5 . 3 , F 5 .0 ) 

READ ( 5, 16,END=90000) ( I LEVEL < N ) * KVI BL < N ) f N=1 ♦ 6 ) 

16 FORMAT ( 6 ( I2»E9.1) ) 

WRITE (6,17) (ILEVEL(N),KVIBL(N),N=1»6) 

17 FORMAT ( ' 1 , 12 »E9.1 ) 

LMAX = 1. + ( ZUPPER - 120/)/DELZ 
READ (5, 18, END =90000) I NCRTM , KPR I M , I PR 

18 FORMAT ( 16, E9. 1,16) 

WRITE (6,19) I NCRTM ,KPRIM, IPR 

19 FORMAT ( ' ' ,I6,1PE9.1,I6) 

READ ( 10,END=90001,ERR=90002) MMM , DELT 1 , LBOUND, DELZ 
X ,DELZ1,LINC, ( Z( ILN ) ,N2Z( ILN ) ,TV( ILN) , ( I I ( I ) , 

X PH I ( I , I LN ) , EPS I L ( I, ILN), NBOLT! I , I LN ) , I = 1 , 6 , 1 ) , I LN = 

X LBOUND, LMAX, LINC ) 

I TEST = MMM 
JMAX = I NCRTM - 1 
DO 1000 M=1 ,10000, INCRTM 
DO 1500 L=LBOUND, LMAX, L INC 
SUM 1 = 0. 

SUM 2 = 0. 

DO 2000 J = 1 ,6 

NLTOT ( J , L ) = NBOLT ( J , L ) + EPSIL(J,L) 

A = (KVIBL(J) - KPRIM)*NBOLT( J,L ) 


00000100 
00000200 
00000300 
00000400 
00000500 
00000600 
00000700 
00000800 
00000900 
00001000 
00001100 
00001200 
00001300 
00001400 
00001500 
00001600 
00001700 
00001800 
00001900 
00002000 
00002100 
00002200 
00002300 
00002400 
00002500 
00002600 
00002700 
00002800 
00002900 
00003000 
0000 3 1 00 
00003200 
00003300 
00003400 
00003500 
00003600 
00003700 
00003800 
00003900 
00004000 
00004100 
00004200 
00004300 
00004400 
00004500 
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B = KVIBL( J)*EPSIL(J,L) 

SUM 1 = SUM1 + A 
2000 SUM2 = SUM 2 + B 

TAUINV = SUM1 + KPRIM*N2Z(L) 

TAUB(L) = 1 ♦ / T AU I NV 
TAUE(L) = 1. /(TAUINV + SUM2) 

TAUN < L > = 1 • / ( 1.3E-12*N2Z(L) ) 

KBQLT(L) = KPRIM + SUM1/N2Z(L) 

KEFF(L) « KBOLT(L) + SUM2/N2Z(L) 

RATIOK(L) = KEFF ( L ) /KBOLT ( L ) 

ALTAUB(L) = AL0G1 0 (TAUB(L) ) 

ALTAUE(L) = AL0G10 ( T AUE ( L ) ) 

ALTAUN(L) = ALOGlO (TAUN { L ) ) 

ALKBOL(L) = ALOGlO ( KBOLT ( L ) ) 

ALKEFF ( L ) = ALOGlO ( KEFF ( L ) ) 

1500 CONTINUE 

: WRITE RESULTS TO TAPE# 

WRITE (11) MMM tDELTl »L BOUND »DELZ » DELZ 1 
X , L I NC , ( Z ( I LN ) , N2Z ( I LN ) * TV( ILN),TAUB( ILN), 

X TAUE ( I LN ) tKBOLT ( I LN ) , KE FF ( I LN ) , R AT I OK < I LN ) , 

X ( 1 1 ( I ) t N60LT ( I f I LN ) t EPS I L ( I v ILN ) •NLTOTU * I LN ) t 

X 1 = 1 ,6, 1 ) 1 1 LN=LBOUNO tLMAX t L INC ) 

IF ( MMM# NE* I TEST ) GO TO 5000 
ITEST = ITEST + IPR 
WRITE ( 6 1 20 ) 

20 FORMAT ( * 1 1 , IX, 1 MMM IS ' ) 

WRITE (6,21) MMM 

21 FORMAT ( • + • , 9X , I 6 ) 

WRITE (6,50) 

50 FORMAT ( 1 0 1 , 2 X , » AL T 1 , 4X , 1 T V 1 ,4X , 1 N ( N2 ) * , 9X , ■ KBOL T 1 9X , 

X 1 KEFF SI IX, ■ RAT I OK 1 ,7X, • TAUN S 10X , • TAUB S 1 1 X , »TAUE» ) 

WRITE ( 6,51 ) ( Z< ILN ) ,TV( ILN) ,N2Z( ILN) , KBOLT ( I LN ) , KEFF (I LN ) , 

X RAT I OK ( I LN ) , T AUN ( ILN),TAUB< ILN),TAUE( I LN ) , I LN = L BOUND , LM AX , L I NC ] 

51 FORMAT ( '0S0PF5.0,0PF6.0,1P7E14#6) 

: PRINTER PLOT. 

XR = 1000. 

XL = 0. 

YT = 11. 

YB = 1. 

NUMPR1 = LMAX + 1 
DO 3000 N = 1 , 6 

DO 3100 LL=LBOUND,LMAX,LINC 
IF (LL.GE.2) GO TO 3105 
ARRAY ( LL ) = 1. 


ARRAY1 CLL ) 


00009000 
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GO TO 3100 

3105 ARRAY ( L L ) = AL0G10< NBOLT ( N * LL ) ) 

ARRAYl(LL) = AL0G10 ( NLTOT ( N * LL ) ) 

3100 CONTINUE 

WRITE (6*22) 

22 FORMAT ( * 1 1 * 40X * 1 DENS I T I E S VERSUS ALTITUDE FOR') 
WRITE (6*23) N 

23 FORMAT ( •+', 71X*I2) 

CALL PL0T2(GRID,XR*XL*YT*YB) 

CALL PL0T3( '*' * Z * ARRAY * NUMPR 1 ) 

CALL PL0T3( ' + ' * Z * ARRAY 1 * NUMP R 1 ) 

CALL PL0T4( 27* 1 LOG NUMBER DENSITY <CM**-3)') 

WRITE (6*24) 

24 FORMAT ( * 0 ' * 59X * ' ALT I TUDE (KM.)*’) 

3000 CONTINUE 

WRITE (6*25) 

25 FORMAT (' 1 • *40X *' REACT ION RATES VERSUS ALTITUDE 1 ) 

YT = -9. 

YB = -14. 

CALL PL0T2(GRID*XR*XL*YT*YB). 

CALL PL0T3( «*' * Z *ALKBOL * NUMPR 1 ) 

CALL PL0T3( • + » * Z * ALKE F * NUMPR 1 ) 

CALL PL0T4(35* 'LOG REACTION RATE ( CM**3 SEC**-1)') 

WRITE (6*24) 

YT = 1.5 
YB = 0.5 
WRITE (6*26) 

26 FORMAT ( ' 1 * * 40X * » R AT I 0 OF REACTION RATES') 

CALL PL0T2 ( GRID*XR*XL*YT *YB) 

CALL PL0T3( * Z * RAT I OK * NUMPR 1 ) 

CALL PL0T4(20* 'RATIO OF KEFF/KBOLT') 

WRITE (6*24) 

YT = 5. 

YB = 0. 

WRITE (6*27) 

27 FORMAT (' 1 '* 40X *' ELECTRON LOSS TIME CONSTANT') 

CALL PLOT 2(GRID,XR*XL*YT *YB) 

CALL PL0T3( '*' * Z * ALT AUB * NUMPR 1 ) 

CALL PL0T3( ' + ' * Z , ALT AUE * NUMPR 1 ) 

CALL PL0T3( ' N ' * Z * ALT AUN * NUM PR 1 ) 

CALL PL0T4124* 'LOG TIME CONSTANT (SEC.)') 

5000 IF ( INCRTM.EO.l ) GO TO 5500 
DO 5500 J J = 1 * JMAX 
READ < 10*END=90001*ERR = 90002 ) 

5500 CONTINUE 
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READ ( 10, END= 90001 , ERR=90002 ) MMM » DELTl , L BOUND, DELZ 
X ,DELZl,UNCt< Z( ILN),N2Z< I LN > , TV ( I LN ) , (II (I) , 

X PHI ( I , I LN ) , EPS I L ( It ILN)tNBOLT( I , I LN ) , I = 1 , 6 , 1 >., I LN = 
X LBOUND , LMAXtLINC ) 

1000 CONTINUE 

END FILE 11 
STOP 

90000 WRITE (6,33) 

33 FORMAT ( 1 1 » , » END OF FILE ENCOUNTERED ON CARD READER') 
STOP 

90001 WRITE (6*84) 

84 FORMAT ( ' 1 • * ' END OF FILE ENCOUNTERED ON INPUT TAPE') 
STOP 

90002 WRITE (6,85) 

85 FORMAT ('l','IO ERROR ENCOUNTERED ON INPUT TAPE') 

STOP 

END 
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